基因家族分析之获取全部cDNA碱基序列构建进化树

生信技能树

发布于 2019-10-09 15:25:40

1.6K0

发布于 2019-10-09 15:25:40

文章被收录于专栏：生信技能树生信技能树

正好，就安排公司学习者以此为课题来一个正规的基因家族分析吧，会有3个教程：

首先是根据基因ID获取基因家族全部cDNA碱基序列构建进化树
然后是根据基因ID获取基因家族全部蛋白质氨基酸序列构建进化树
最后是根据基因家族的motif查询全部基因后获取序列构建进化树

现在让我们看第一个教程，以趋化因子基因家族为背景来获取序列进行多序列比对后绘制系统发育树。

PS：本教程绝大部分都是公司学习者辛苦摸索写的，只有一个perl单行命令是我操刀。

趋化因子背景知识

趋化因子被分为四个主要的亚家族：CXC，CC，CX3C和XC。所有这些蛋白质都通过与G蛋白质连接的跨膜受体（称为趋化因子受体）相互作用而发挥其生物学作用，该受体在其靶细胞表面被选择性地发现。

维基百科关于趋化因子介绍：https://zh.wikipedia.org/wiki/

参考文章：Charo IF, Ransohoff RM. The many roles of chemokines and chemokine receptors in inflammation. N Engl J Med. 2006 Feb 9;354(6):610-21.

也称做趋化激素、趋化素或是化学激素。

是一小分子细胞因子家族蛋白。趋化因子蛋白的共同结构特征包括，分子量小（约8-10 千道尔顿），有四个位置保守的半胱氨酸残基以保证其三级结构。这些小蛋白因其有定向细胞趋化作用而得名。当然，这些蛋白有些趋化因子历史上还有其他的名字，包括已知的SIS细胞因子家族、 SIG细胞因子家族，SYC细胞因子家族和血小板因子-4家族。有的趋化因子被认为促进炎症反应，而有些趋化因子被认为在正常的修复过程或发育中控制细胞的迁徙。在所有脊椎动物和一些病毒和一些细菌中有趋化因子存在，但不存在于其他无脊椎动物。这些蛋白质结合到趋化因子受体而起作用，趋化因子受体是G蛋白偶联受体，选择性地表达在靶细胞表面。

当然，命名往往不是那么的规律，比如IL-8是其中一个成员，在该基因家族中的ID是CXCL8。

典型结构：

数据下载

首先要拿到CCL基因家族的基因列表，可以在人类基因组命名委员会数据库下载： https://www.genenames.org/data/genegroup/#!/group/483
去gencode拿到数据库下载最新的gtf文件 https://www.gencodegenes.org/human/release_30.html wget -c ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/gencode.v32.annotation.gtf.gz
去ensemble数据库下载人类的cdna的序列 wget -c ftp://ftp.ensembl.org/pub/release-97/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz

家族成员一览表

前面我们介绍过趋化因子被分为四个主要的亚家族：CXC，CC，CX3C和XC，全部基因下载如下：

#### CC chemokines

Name    Gene    Other name(s)   Receptor    Uniprot
CCL1    Scya1   I-309, TCA-3    CCR8    
CCL2    Scya2   MCP-1   CCR2    P13500
CCL3    Scya3   MIP-1a  CCR1    P10147
CCL4    Scya4   MIP-1β  CCR1, CCR5  P13236
CCL5    Scya5   RANTES  CCR5    P13501
CCL6    Scya6   C10, MRP-2  CCR1    P27784
CCL7    Scya7   MARC, MCP-3 CCR2    P80098
CCL8    Scya8   MCP-2   CCR1, CCR2B, CCR5   P80075
CCL9/CCL10    Scya9   MRP-2, CCF18, MIP-1?    CCR1    P51670
CCL11    Scya11  Eotaxin CCR2, CCR3, CCR5    P51671
CCL12    Scya12  MCP-5       Q62401
CCL13    Scya13  MCP-4, NCC-1, Ckβ10 CCR2, CCR3, CCR5    Q99616
CCL14    Scya14  HCC-1, MCIF, Ckβ1, NCC-2, CCL   CCR1    Q16627
CCL15    Scya15  Leukotactin-1, MIP-5, HCC-2, NCC-3  CCR1, CCR3  Q16663
CCL16    Scya16  LEC, NCC-4, LMC, Ckβ12  CCR1, CCR2, CCR5, CCR8  O15467
CCL17    Scya17  TARC, dendrokine, ABCD-2    CCR4    Q92583
CCL18    Scya18  PARC, DC-CK1, AMAC-1, Ckβ7, MIP-4       P55774
CCL19    Scya19  ELC, Exodus-3, Ckβ11    CCR7    Q99731
CCL20    Scya20  LARC, Exodus-1, Ckβ4    CCR6    P78556
CCL21    Scya21  SLC, 6Ckine, Exodus-2, Ckβ9, TCA-4  CCR7    O00585
CCL22    Scya22  MDC, DC/β-CK    CCR4    O00626
CCL23    Scya23  MPIF-1, Ckβ8, MIP-3, MPIF-1 CCR1    P55773
CCL24    Scya24  Eotaxin-2, MPIF-2, Ckβ6 CCR3    O00175
CCL25    Scya25  TECK, Ckβ15 CCR9    O15444
CCL26    Scya26  Eotaxin-3, MIP-4a, IMAC, TSC-1  CCR3    Q9Y258
CCL27    Scya27  CTACK, ILC, Eskine, PESKY, skinkine CCR10   Q9Y4X3
CCL28    Scya28  MEC CCR3, CCR10 Q9NRJ3

#### CXC chemokines

Name    Gene    Other name(s)   Receptor    Uniprot
CXCL1    Scyb1   Gro-a, GRO1, NAP-3, KC  CXCR2   P09341
CXCL2    Scyb2   Gro-β, GRO2, MIP-2a CXCR2   P19875
CXCL3    Scyb3   Gro-?, GRO3, MIP-2β CXCR2   P19876
CXCL4    Scyb4   PF-4    CXCR3B  P02776
CXCL5    Scyb5   ENA-78  CXCR2   P42830
CXCL6    Scyb6   GCP-2   CXCR1, CXCR2    P80162
CXCL7    Scyb7   NAP-2, CTAPIII, β-Ta, PEP       P02775
CXCL8    Scyb8   IL-8, NAP-1, MDNCF, GCP-1   CXCR1, CXCR2    P10145
CXCL9    Scyb9   MIG, CRG-10 CXCR3   Q07325
CXCL10    Scyb10  IP-10, CRG-2    CXCR3   P02778
CXCL11    Scyb11  I-TAC, β-R1, IP-9   CXCR3, CXCR7    O14625
CXCL12    Scyb12  SDF-1, PBSF CXCR4, CXCR7    P48061
CXCL13    Scyb13  BCA-1, BLC  CXCR5   O43927
CXCL14    Scyb14  BRAK, bolekine      O95715
CXCL15    Scyb15  Lungkine, WECHE     Q9WVL7
CXCL16    Scyb16  SRPSOX  CXCR6   Q9H2A7
CXCL17    VCC-1   DMC, VCC-1      Q6UXB2

#### C chemokines

Name    Gene    Other name(s)   Receptor    Uniprot
XCL1    Scyc1   Lymphotactin a, SCM-1a, ATAC    XCR1    P47992
XCL2    Scyc2   Lymphotactin β, SCM-1β  XCR1    Q9UBD3

#### CX3C chemokines

Name    Gene    Other name(s)   Receptor    Uniprot
CX3CL1    Scyd1   Fractalkine, Neurotactin, ABCD-3    CX3CR1  P78423

有了如上所述的基因信息，就可以根据我们从gencode下载的gtf文件进行。

获取CCL家族基因的转录本id

在人类基因组命名委员会数据库下载拿到的是一个csv文件，打开文件，第二列Approved symbol就是gene symbol id，总共有45个，保存到一个文件中，如：CCL_genes

$ cat CCL_genes
CCL1
CCL2
CCL3
CCL3L1
CCL3L3
......

然后从gtf文件中拿到CCL基因家族的所有转录本id，共124个，以及对应的蛋白质id（部分转录本没有蛋白id）

$ zcat gencode.v32.annotation.gtf.gz |grep -f CCL_genes|awk 'BEGIN{FS="\t";OFS="\t"}{if($3=="transcript")print $0}'|cut -d ";" -f 2,4,6,8|less >CCL_tmp.txt

# 可以进一步做以下处理，比较复杂，总之就是拿到转录本id，基因名，和对应的蛋白质id
$ less CCL_tmp.txt |sed 's/transcript_support_level..../ /' | cut -d " " -f 3,5,7,9 | sed 's/\"//g' |tr ';' '\t' >CCL_trans_gene_pro_id.txt

$ less CCL_trans_gene_pro_id.txt
ENST00000367819.3        XCL2    XCL2-201        ENSP00000356793.2
ENST00000367818.4        XCL1    XCL1-201        ENSP00000356792.3
ENST00000646475.1        CCL20   CCL20-205       ENSP00000496658.1
ENST00000358813.5        CCL20   CCL20-201       ENSP00000351671.4
ENST00000409189.7        CCL20   CCL20-202       ENSP00000386273.3

（PS：我觉得上面的grep脚本运行肯定会很慢，但是我懒得去修改它了。）

可以把转录本id单独保存为一个文件，比如我这里保存为CCL_trans_id.txt

$ cat CCL_trans_id.txt
ENST00000367819.3
ENST00000367818.4
ENST00000646475.1
ENST00000358813.5
ENST00000409189.7
ENST00000489160.1
......

获取转录本的碱基序列

从人类cdna的fasta文件中提取提取上面转录本id的碱基序列这里需要用到一个perl 单行命令：（只有这个命令是我写的）

cat ccl_t.id Homo_sapiens.GRCh38.cdna.all.fa |perl -alne '{if(/^>/){if(exists($h{$F[0]})){$tmp=1}else{$tmp=0}};print if $tmp==1;$h{">$_"}=1}'

使用起来其实很简单，需要两个输入文件，第一个是要获取的转录本id号，即上面的CCL_trans_id.txt，第二个是人类的Homo_sapiens.GRCh38.cdna.all.fa，结果输出就是fasta格式，可以保存为CCL_trans_cdna.fa

转录本id转成CCL

把id替换成CCL，这样子可读性比较强一些，处理起来有点复杂，首先我们先把上面拿到的CCL_trans_cdna.fa，每一条序列合为一行，然后用awk命令将处理后的结果与前面的CCL_trans_gene_pro_id.txt按指定行进行合并，再进行一定的处理，最后转成fasta格式，比较复杂，需要花较多的时间理解下面的代码：

$ less CCL_trans_cdna.fa |tr '\n' '\t'|tr '>' '\n' >1.tmp

$ cat CCL_trans_gene_pro_id.txt 1.tmp |awk '{a[$1]=a[$1]"\t"$0}END{for( i in a){print ">"a[i]}}'|cut -f 1,4,7-|sed 's/ //' |sed 's/\t//'|awk '{if($2!="")print $0 }'|tr '\t' '\n'|sed '/^$/d' >all_id.fa

做这么多操作就是为了一个目的：将CCL_trans_cdna.fa 中的id转换成可读性较强的基因id

$ less CCL_trans_cdna.fa |grep '>'
>ENST00000557161.1
>ENST00000259607.7
>ENST00000378792.1
>ENST00000311925.7
.....

$ less all_id.fa |grep '>'
>CCL4L2-210
>CCL4-203
>CCL22-201
>CCL3L1-202
......

拿到最后的all_id.fa文件后，就可以进行进化树分析了

绘制CCL基因家族的cdna进化树

library(Biostrings)

## merge
## fasta import & auto-masking
merge_8 <- readAAStringSet("all_id.fa", 
                           format="fasta",
                           nrec=-1L, 
                           skip=0L, 
                           seek.first.rec=FALSE, 
                           use.names=TRUE)
aln <- muscle::muscle(merge_8)
auto <- maskGaps(aln, min.fraction=0.5, min.block.width=4)
## get the tree!
sdist <- stringDist(as(auto,"AAStringSet"), method="hamming")
clust <- hclust(sdist, method = "single")
png(paste0('HM_CCL_tree', '.png'),res = 150,width = 1560,height = 1080)
plot(as.phylo(clust),type="fan",cex = 1.0,main="HM CCL Gene Group Pylogenetic Tree")
dev.off()