每日学习20200922

Posted on

基因的跨物种序列比对以及进化树构建

工具：muscle，ggtree

1. 序列比对工具

Muscle：本地版和网页版都有，

本地版：https://www.drive5.com/muscle/downloads.htm

网页版： https://www.ebi.ac.uk/Tools/msa/muscle/

参考网站： https://cloud.tencent.com/developer/article/1625177

ClustalW：速度较慢

MAFFT：氨基酸的多序列比对

lastZ:基因组比对，适用于染色质序列比对，对于gap友好

T-coffee：整合多种信息，比如结构和实验，功能比较强大

PHAST: 自荐

2. 准备序列

在ncbi下载数据，注意是基因组还是转录组的fastq数据。

如果要导入snapview提前查看序列及其特征，需要额外下载ncbi genebank 格式的文件。

各物种的序列将其整合到一个文件中，需要注意抬头的书写，识别的时候会根据空格进行分割，所以特别注意序列名称的特异性。

3. 过程：

    1. 首先，使用网页版的muscle
      INPUT: 上述的fasta序列，在output format参数选择中，选择fasta格式和clustalW格式
      OUTPUT：jalview launch file可以导入到ebi推荐的桌面app中查看 ，alignment格式文件另存为，简单的进化树直接生成。
      1. 在R中绘制进化树，与序列对齐，使用ggtree。

      特别注意

      muscle默认输出的ph格式是ntk格式，可以用ggtree的read.newick 读取。

      fasta的抬头需要和ph文件的抬头一致。

    1. 如果进化树的序列知显示一部分，则需要对fasta输出文件进行截取。

代码参考：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
library(ggtree)
library(treeio)
setwd("/media/ggj/NEW/DifferentiationForce/xbp1/xbp1-sequnece/species/muscle-mRNA-result/")

## from muscle website file
tree<-read.newick("muscle-I20200922-083912-0228-12080513-p1m.ph")
tree

## from muscle website file
fasta<-"/media/ggj/NEW/DifferentiationForce/xbp1/xbp1-sequnece/species/muscle-mRNA-result/aln-fasta-muscle-I20200922-074229-0213-98387757-p1m_rename.fasta"
x<-read.fasta(fasta)

## check the taxa names
setdiff(tree$tip.label,labels(x))
setdiff(labels(x),tree$tip.label)


#p3 <- ggtree(tree,) +geom_tiplab(hjust = -0.05,size=6,fontface="italic") + xlim(NA,1.2)

p3 <- ggtree(tree)+geom_tiplab(hjust = -0.05,size=3,fontface="italic")
print(p3)
fasta<-"/media/ggj/NEW/DifferentiationForce/xbp1/xbp1-sequnece/species/muscle-mRNA-result/aln-fasta-muscle-I20200922-074229-0213-98387757-p1m_rename.fasta"
library(RColorBrewer)
col<-brewer.pal(5,"Paired")
msaplot(p3,fasta = fasta,offset=0.05,color = col)



############# subset nls and lecine zipper sequence
## mannal select the 600-900bp region
fasta<-"aln-fasta-muscle-I20200922-074229-0213-98387757-p1m_rename_part.fasta"
xx<-read.fasta(fasta)

p3 <- ggtree(tree)+geom_tiplab(hjust = -0.05,size=3,fontface="italic")
#print(p3)
library(RColorBrewer)
col<-brewer.pal(5,"Paired")
msaplot(p3,fasta = xx,offset=0.05,color = col)

20200815 每日学习

文件格式： bigWig, bed, bigBed, bed narrowPeak

http://genome.ucsc.edu/goldenPath/help/wiggle.html

bigWig: