研究背景

随着高通量测序技术的迅猛发展，生命树各主要进化分支组学数据不断完善，使得从分子水平解析宏观进化生物学问题成为可能。通过转录组测序获得的基因表达变化的差异往往蕴含着表型进化的信息，因而追踪并量化转录组在进化过程中的变化规律，将为我们理解转录组变化如何驱动表型进化提供新的认知 （Domazet-Lošo and Tautz, 2010）。Domazet-Lošo等（2007）提出了系统进化层级理论，认为基因组中每个基因都有其系统进化层级，并可通过与已知基因组进行相似性比对来确定基因的“进化年龄”。基于这一理论所建立的转录组年龄指数（TAI）分析方法，采用基因表达量的基因年龄加权总和除以基因表达量的总和得到的TAI指数来表征转录组年龄，即TAI值越高代表转录组越年轻 （Domazet-Lošo and Tautz, 2010）。该方法已广泛应用于动物、植物、微生物等类群，在发育沙漏模型进化推演、干细胞起源进化、海洋幼虫起源进化等重要科学问题解析中取得突破性发现（Domazet-Lošo and Tautz, 2010; Hemmrich, et al., 2012; Futo et al., 2020; Wang et al., 2020）。在Domazet-Lošo等的研究基础上，Hajk-Georg等（2018）开发了可用于TAI分析的R包，进一步为研究者提供了一个标准化、自动化的TAI分析流程。但现有的TAI分析方法在解析宏观进化生物学问题时仍存有局限性，如Domazet-Lošo等开发的算法仅计算各系统进化层级对进化终端（即现存物种）转录组TAI的贡献，不能计算各祖先状态转录组的TAI，故无法实现对祖先转录组变化规律的推演分析。针对已有TAI算法的分析局限性，本方案在Domazet-Lošo等开发的原始算法的基础上，提出了一种改进的TAI分析方法（Wang et al., 2020），通过逐层移除相对更年轻的基因来推演各祖先进化节点的转录组状态，并利用该方法实现了对转录组宏观进化规律的推演。

应用程序和软件

1. Perl (https://www.perl.org/)
2. fqtrim (https://ccb.jhu.edu/software/fqtrim/)
3. STAR (https://github.com/alexdobin/STAR)
4. Samtools (http://www.htslib.org/)
5. featureCounts (http://bioinf.wehi.edu.au/featureCounts/)
6. Trinity (https://github.com/trinityrnaseq/trinityrnaseq)
7. TransDecoder (https://github.com/TransDecoder/TransDecoder)
8. BLAST+ (ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/ )
9. R (https://www.r-project.org/)
10. myTAI (https://github.com/drostlab/myTAI)
    

数据分析

一、基因表达谱计算

1. 从NCBI网站的SRA数据库(https://www.ncbi.nlm.nih.gov/sra/)中下载所需的高通量测序数据（表1）。
   
   表1 虾夷扇贝不同发育时期的转录组数据
   
   
   
2. 测序数据预处理
   对高通量测序平台测序获得的原始数据进行质量过滤，用fqtrim软件去掉含不确定碱基N的序列以及低质量reads（超过20%的碱基质量值小于20的reads）,获得的clean reads用于后续分析，输入文件为FASTQ格式。
   $ fqtrim -A -q 20 -w 80 -m 0 -p 8 -o trimed.fastq sample.R1.fastq, sample.R2.fastq -P 33
   
3. 表达谱计算（无参照基因组）
   3.1

   对于没有参照基因组的物种，首先需要用Trinity软件将预处理获得的clean reads进行de novo拼接。
   $ Trinity.pl --seqType fq –max_memory 40G --left sample1.R1.clean.fq.gz, sample2.R1.clean.fq.gz –right sample1.R2.clean.fq.gz, sample2.R2.clean.fq.gz --CPU 16 --full_cleanup > test.log
   注：--left : left reads，多个文件以逗号隔开
           --right : right reads，多个文件以逗号隔开
          Trinity运行过程中会自动检测之前输入，如果中途发生报错，进行修改之后，会在原有输出文件的基础上继续运行，不会对之前结果进行覆盖
   
   3.2

   对拼接获得的isoform聚类，每条unigene保留最长的转录本，以此作为转录组比对的参照数据库
   $ perl trinity_unigene.pl -i trinity.fasta
   
   3.3

   用TransDecoder软件对3.2获得的数据库进行结构注释
   $ TransDecoder.LongOrfs -t trinity_unigene.fasta
   $ TransDecoder.Predict -t trinity_unigene.fasta --single_best_only
   
4. 表达谱计算（有参照基因组）
   4.1

   用STAR软件建立参照基因组的索引文件，将结果文件夹命名为genome_index（图1）。
   4.2

   $ STAR --runThreadN 8 --runMode genomeGenerate –genomeDir ./genome_index/ --genomeFastaFiles genome.fasta --sjdbOverhang 149 --genomeChrBinNbits 14 > starindex.log
   注：--genomeDir：需要预先建立一个有可读写权限的文件夹，用于存放生成的索引文件
          --sjdbOverhang：该数值设定为单端reads的长度减1，是在注释可变剪切序列时使用的最大长度值
          --genomeChrBinNbits：默认值为18，如需节约运行内存，可根据参考基因组的大小重新定义，=min(18, log2(GenomeLength/NumberOfReference
   s))
   
   
   图1 genome_index文件夹中的索引文件
   
   
   4.3

   将预处理获得的clean reads用STAR软件进行比对，将结果文件夹命名为STAR_out（图2）
   $ STAR --runThreadN 8 –genomeDir ./genome_index / --readFilesCommand zcat --readFilesIn R1.clean.fastq.gz R2.clean.fastq.gz --sjdbGTFtagExonParentGene genome.gff3 --outFilterMultimapNmax 100000 --outFileNamePrefix ./STAR_out/ --outReadsUnmapped Fasta > star.log
   注：--outFilterMultimapNmax：默认值为100000，限制SAM文件的大小，可根据数据量大小适当调高该参数。
          --FileNamePrefix：需要事先建立一个有可读写权限的文件夹，用于存放输出文件
   
   
   图2. STAR_out文件夹中的结果文件
   
   
   4.4

   将结果文件夹下的sam文件转换为bam文件
   $ samtools view -bS ./STAR_out/Aligned.out.sam > sample.bam
   
   4.5

   根据基因组结构注释文件，用featureCounts软件统计map到各基因外显子区的reads数
   $ featureCounts -p -t exon -g Parent -a genome.gff3 -o sample.txt sample.bam
   注：-p：只能用在双端测序paired-end的情况中，会统计fragment而不统计read
          -t：根据gff文件设置feature-type，同时read只有落到这些feature上才会被统计到，默认是“exon”。
          -g：根据gff文件设置id identifier，要确保-g 参数指定的值是gff文件中的唯一标识符，默认是“gene_id”。
   
   4.6

   将基因的count值转换为TPM 值
   $ perl count2TPM.pl -l exon-length -i count.txt -o TPM.txt
   注：-l：步骤3.4中定义的feature的长度，输入文件包括两列，第一列为基因ID（geneID），第二列为feature长度（geneL），两列数据以tab键分隔。
          -i：步骤3.4中统计获得的count值，输入文件第一列为基因ID（geneID），第二列及之后为count值（sample1, sample2…），每列数据以tab键分隔。
   

二、定义基因年龄

基因年龄的定义方法参照Domazet-Lošo et al. (2010)开发的系统层级法(Phylostratigraphy)。

1. 定义系统进化层级
   根据NCBI网站Taxonomy数据库（https://www.ncbi.nlm.nih.gov/taxonomy）中的种系发生树，将全基因组的基因划分为不同的系统进化层级（图3，图4）。以虾夷扇贝为例，本研究中将虾夷扇贝的基因划分为11个进化层级（图5），其中将细胞生物（cellular organisms）定义为系统进化层级的第一级（PS1），真核生物（Eukaryota）定义为系统进化层级的第二级（PS2），依次类推，将虾夷扇贝属（Mizuhopecten）定义为系统进化层级的最后一级（PS11）。
   
   
   图3. 虾夷扇贝在Taxonomy数据库中的进化层级
   
   
   图4. NCBI的Protein数据库中属于Mollusca进化层级的蛋白序列
   
   
   图5. 本研究中将虾夷基因划分为11个进化层级（PS1-PS11）
   
   
2. 定义基因年龄
   将虾夷扇贝基因组的全部基因与NCBI的Protein数据库做同源比对，以基因能比对到的最古老的物种所在的进化层级定义基因年龄，在任何物种中都比不到的定义为虾夷扇贝特有的新基因，划为PS11。同源比对的命令如下：
   $ blastp -db database.fasta -query genome_protein.fasta -evalue 1e-3 -out output.txt -max_target_seqs 1 -num_threads 8 -outfmt 6
   

三、转录组年龄指数（TAI）分析

参考Wang et al., (2020)发表的文章中关于TAI分析的计算公式，使用R（version 3.5.0）软件中的myTAI 软件包进行TAI相关分析。

1. 通过CRAN镜像安装myTAI (0.8.0)软件包
   >source("http://bioconductor.org/biocLite.R")
   >biocLite('myTAI')
   
2. 整理TAI分析所需的输入文件
   输入文件的第一列是基因对应的进化层级，第二列是GeneID，第三列之后是基因在各发育时期的表达量，如下图所示。输入文件保存为csv格式（图6）。
   
   
   图6. TAI输入文件格式
   
   
3. 调用myTAI 软件包进行TAI计算（图7）
   >library(myTAI)
   >data<-read.csv('TAI输入文件.csv', header=T)
   >dataset<-as.data.frame(data)
   >TAI(dataset)
   
   
   图7. 虾夷扇贝16个发育时期的global TAI值
   
   注：这里需要着重强调的是，根据TAI计算公式，
   
   ps_i 表示基因i的进化层级，e_is 表示基因i在发育时期s时的表达量，n 表示用于TAI计算的基因总数。
   当计算虾夷扇贝某一发育阶段的global TAI时，输入文件中包含PS1-PS11的全部基因。在计算祖先时期TAI值时，我们的计算方法与Domazet-Lošo and Tautz (Domazet-Lošo and Tautz, 2010)不同。我们通过逐层去除相对更年轻的基因，来推演祖先时期的转录组状态。例如，要计算扇贝转录组在双壳纲祖先 (PS10) 时的TAI，则在输入文件中完全去除PS11基因。以此类推，在计算处于进化层级N时的TAI值时，起源晚于N的基因总是从输入文件中移除（即，从计算公式的“分子”和“分母”中同时移除）。
   
4. 发育时期TAI变化曲线及显著性分析（图8，图9）
   >PlotSignature (ExpressionSet=dataset,
                  measure="TAI",
                  TestStatistic="FlatLineTest",
                  xlab="Ontogeny",
                  ylab="TAI")
   用于显著性分析的FlatLineTest功能其具体参数如下：
   FlatLineTest (ExpressionSet = dataset,
                 permutations = 1000,
                 plotHistogram = TRUE )
   
   
   图8. 虾夷扇贝发育时期global TAI变化曲线及显著性分析
   
   
   图9. 虾夷扇贝处于不同进化层级时的TAI曲线
   
   
5. 计算每个基因对global TAI的贡献值，输出csv格式的结果文件（图10）
   >partialTAI<-pMatrix(dataset)
   >head(partialTAI)
   >write.csv(partialTAI,'pMatrix.csv')
   注：根据TAI计算公式，所有基因在某一发育时期的TAI贡献值之和等于该时期的global TAI值。
   
   
   图10. 虾夷扇贝各基因对global TAI的贡献值
   
   
6. 计算不同进化层级对TAI的相对贡献（图11）
   进化层级N中的基因对转录组TAI的相对贡献用PS_S(N) 表示，计算公式为
   
   这里，表示属于进化层级N的所有基因在发育时期s中的表达量之和，e_is表示基因i在发育时期s时的表达量，n 表示进化层级N的基因数。
   
   
   图11. 不同进化层级基因对扇贝担轮幼虫TAI的相对贡献
   
   

致谢

本研究受国家自然科学基金项目(31871499 & 31900369)，山东省泰山学者计划项目和山东省博士后创新项目资助。

参考文献

1. Domazet-Lošo, T., Brajković, J. and Tautz, D. (2007). A phylostratigraphy approach to uncover the genomic history of major adaptations in metazoan lineages. Trends in Genetics, 23(11):533-539.
2. Domazet-Lošo, T. and Tautz, D. (2010). A phylogenetically based transcriptome age index mirrors ontogenetic divergence patterns. Nature 468(7325): 815-818.
3. Drost, H.G., Gabel, A., Liu, J., Quint, M. and Grosse, I. (2018). myTAI: evolutionary transcriptomics with R. Bioinformatics, 34(9):1589-1590.
4. Futo, M., Opašić, L., Koska, S., Čorak, N., Široki, T., Ravikumar, V., Thorsell, A., Lenuzzi, M., Kifer, D., Domazet-Lošo, M., Vlahoviček, K., Mijakovic, I. and Domazet- Lošo, T. (2020). Embryo-like features in developing Bacillus subtilis biofilms. Molecular Biology and Evolution, 38(1):31-47.
5. Hemmrich, G., Khalturin, K., Boehm, A.M., Puchert, M., Anton-Erxleben, F., Wittlieb, J., Klostermeier, U.C., Rosenstiel, P., Oberg, H.H., Domazet-Lošo, T., Sugimoto, T., Niwa, H. and Bosch, T.C.G. (2012). Molecular signatures of the three stem cell lineages in Hydra and the emergence of stem cell function at the base of multicellularity. Molecular Biology and Evolution, 29(11):3267-3280.
6. Wang, J., Zhang, L., Lian, S., Qin, Z., Zhu, X., Dai, X., Huang, Z., Ke, C., Zhou, Z., Wei, J., Liu, P., Hu, N., Zeng, Q. and Dong, B. (2020). Evolutionary transcriptomics of metazoan biphasic life cycle supports a single intercalation origin of metazoan larvae. Nature Ecology & Evolution, 4(5): 725-736.