土壤宏转录组学样本前处理与数据分析
Sample Pretreatment and Data Analysis of Soil Metatranscriptome   

材料与试剂

1. 50 ml离心管 (Thermo Fisher Scientific, USA)
2. 酚氯仿异戊醇 (配比25：24：1, pH=8)
3. 琼脂糖
4. 电泳缓冲液 (TAE)
5. RNA提取试剂盒RNA Mini Kit (Qiagen, Hilden, Germany)
   

仪器设备

1. 超微量紫外分光光度计NanoDrop One (Thermo Fisher Scientific, MA, USA)
2. 台式高速冷冻离心机 (Techcomp (Holdings), model: CT15RT)
3. 涡旋仪 (Qiangen, catalog number: 13000-V1-15)
4. 水浴锅
5. 移液枪
6. 电泳仪
   

数据分析软件及平台

1. CD-HIT (http://weizhongli-lab.org/cd-hit/)
2. Fastp (v0.19.5) (https://github.com/OpenGene/fastp)
3. SortMeRNA (v2.1) (https://github.com/biocore/sortmerna)
4. IDBA_tran (v1.1.3) (https://github.com/loneknightpy/idba)
5. Prodigal (v.2.6.3) (https://github.com/hyattpd/Prodigal)
6. CD-HIT (v4.8.1) (http://weizhongli-lab.org/cd-hit/)
7. RSEM (v1.3.1) (https://github.com/deweylab/RSEM)
8. Bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)
9. Diamond (v0.1.10) (http://www.diamondsearch.org)
   注：上述软件均在Linux操作系统下进行。
   

实验步骤

一、取样过程
采集0-20 cm湿地土壤样品5 g于50 ml的无菌离心管中，加入RNA抑制酶抑制剂约10 ml使其浸没土壤样品并封存。低温运输并尽快于-80 °C保存。

二、注意事项
确保实验工作区无RNase污染并且整个操作过程戴橡胶手套。

三、RNA制备

1. 加2 g土壤到15 ml磁珠管 (试剂盒提供) 中。
2. 依次向离心管内加入2.5 ml Bead solution溶液，0.25 ml SR1溶液和0.8 ml的IRS溶液并漩涡混匀。
   发生的反应：Bead Solution是一种缓冲液可用来打散细胞和土壤颗粒；SR1能帮助细胞裂解，可以破坏脂肪酸和几种微生物细胞膜相关的脂类；IRS可以帮助去除腐殖质、细胞碎片和蛋白质等杂质。
3. 向磁珠管加入3.5 ml酚氯仿异戊醇溶液 (pH = 8)，漩涡混匀直到分层消失。
4. 最大转速涡旋混匀15 min。
   发生反应：从1到4步的化学试剂和漩涡使细胞裂解，酚/氯仿/异戊醇使其最大程度裂解，溶解的细胞和试剂混合在一起，蛋白质降解只剩下核酸在溶液中。
5. 2,500 × g 离心10 min。
   发生反应：离心导致混合样品相分离。离心后能观察到三相，底下的有机相包括蛋白质和细胞碎片，中间相包括腐殖质和其他有机及无机物质，上层相包括所有的核酸。
6. 小心转移上层水相于一新15 ml离心管中 (试剂盒提供)。
   注：用枪头小心吸取上层水相，误触碰界面。
7. 加1.5 ml SR3到水相中，漩涡混合。4 °C孵育10 min， 2,500 × g离心10 min。
8. 将上清液转移到一新15 ml离心管中 (不要碰到下面的沉淀物)，加入5 ml SR4混匀，室温放置30 min。
   注：分层明显的界面处小心吸取上清，勿刺破 (触碰) 界面。
9. 2,500 × g离心30 min。
10. 倒出上清液，将离心管倒置在纸巾上5 min。
    注：依据土壤类型的不同，沉淀可能较大或颜色较深 (Mettel et al., 2010)。
11. 摇晃SR5溶液使其混合，加1 ml SR5溶液到离心管中，使沉淀再完全悬浮。
    注：沉淀可能由于土壤样品的不同不易悬浮，可能需要将离心管放到45 °C的水浴池中10 min再悬浮，再漩涡混合，重复这样直到沉淀重悬浮。
12. 为每个RNA样品准备一个捕集柱。
    12.1

    将捕集柱悬挂到离心管上。
    12.2

    加2 ml SR5溶液到捕集柱上，使其重力流。允许SR5溶液完全流过捕集柱。
13. 将11步的RNA分离样加到捕集柱中，使其流过捕集柱。
14. 用1 ml SR5溶液洗涤捕集柱，流出液收集在15 ml离心管中。
    反应：样品加到捕集柱上，核酸结合到柱基质上。捕集柱用SR5溶液洗涤确保未结合的污染物被去除掉。
15. 将捕集柱转移到一新15 ml离心管中，摇晃SR6，然后加1 ml SR6溶液到捕集柱中使其流过捕集柱，洗提RNA。
    发生反应：SR6溶液RNA洗提缓冲液是专有的盐溶液，它能使RNA流出而DNA、剩下的细胞碎片和抑制剂依然留在捕集柱上。
16. 将洗提的RNA转移到2.2 ml离心管中，并加1 ml SR4，至少倒置混合一次，-20 °C静置10 min。
17. 13,000 × g离心15 min。
18. 移去上清液，将RNA离心管倒置在纸巾上10 min，风干颗粒物。
19. 加100 μl SR7溶液使RNA颗粒再悬浮。
    

结果与分析

一、土壤RNA样本提取效果评价
实验借助市售土壤RNA提取试剂盒 (RNA Mini Kit, Qiagen, Germany) 比较两种不同类型 (酸性、碱性) 的湿地土壤样本总RNA提取效果，纯度和浓度测试结果如表1所示。OD260代表核酸的吸光度，OD280代表蛋白质的吸光度，OD230代表其他杂质 (多糖等) 的吸光度。OD是光密度值。一般来说，OD260/OD280介于1.9~2.0 说明RNA提取纯度高，污染小。OD260/OD280 < 1.7时表明有蛋白质或酚污染；OD260/OD280 > 2.0时表明可能有异硫氰酸残存。原核生物的核糖体rRNA主要由23S、16S和5S rRNA组成。实验室主要通过观察核糖体的23S rRNA和16S rRNA条带的亮度和片段形状来判定RNA的提取效果 (Peano et al., 2013)，本实验中根据OD260/OD280介于1.9~2.0之间，16S rRNA 及 23S rRNA两条标志性条带清晰，说明该方法提取RNA得到了比较纯的RNA样品，经公司评价，满足建库要求。将装有RNA样品的EP管用干冰密封好，送至美格基因进行宏转录组测序 (测序平台为Illumina Hiseq Xten)。

表1. 部分样本RNA提取浓度和纯度

图1. 两种类型湿地土壤样品试剂盒提取的总RNA琼脂糖凝胶电泳图。1-4分别代表酸性湿地土壤RNA，5-8分别代表碱性湿地土壤RNA，M代表片段大小不同的核酸标记物。

二、宏转录组下机数据分析流程

1. 测序得到的原始数据 (raw data) 中包含接头序列及低质量碱基 (Q < 30)，首先经过fastp v0.19.5 (https://github.com/OpenGene/fastp) 去接头 (adapter) 及过滤碱基质量后得到高质量序列 (clean data) 以便用于后续分析：
   # Installation
   $ conda install -c bioconda fastp
   # Quality filter
   $ fastp -I $Sample_R1.fq.gz -I $Sample_R2.fq.gz -o $Sample_clean_R1.fq -O $Sample_clean_R2.fq -l 50 -q 30 -t 10
2. 高质量序列中仍包含大量的核糖体RNA (rRNA)，通过SortMeRNA version 2.1b (https://github.com/biocore/sortmerna) 过滤clean data中的rRNA序列：
   # Installation
   $ conda install -c bioconda sortmerna=2.1b
   # rRNA filter
   $ sortmerna --ref $refdir/rRNA_databases/silva-arc-16s-id95.fasta, $refdir /index/silva-arc-16s-id95: \
   $refdir/rRNA_db/silva-arc-23s-id98.fasta, $refdir/index/silva-arc-23s-id98: \
   $refdir/rRNA_db/silva-bac-16s-id90.fasta, $refdir/index/silva-bac-16s-id90: \
   $refdir/rRNA_db/silva-bac-23s-id98.fasta, $refdir/index/silva-bac-23s-id98: \
   $refdir/rRNA_db/silva-euk-18s-id95.fasta, $refdir/index/silva-euk-18s-id95: \
   $refdir/rRNA_db/silva-euk-28s-id98.fasta, $refdir/index/silva-euk-28s-id98: \
   $refdir/rRNA_db/rfam-5s-database-id98.fasta, $refd-ir/index/rfam-5s-database-id98: \
   $refdir/rRNA_db/rfam-5.8s-database-id98.fasta, $refd-ir/index/rfam-5.8s-database-id98 \
   --reads $Sample_clean_R1.fq --reads $Sample_clean_R2.fq --fastx --other $Sample_rmRNA --aligned $Sample_aligned -a 30 -paired-out
   $ unmerge-paired-reads.sh $Sample_rmRNA.fq $Sample_rmRNA_R1.fq $Sample_rmRNA_R2.fq
   rRNA数据库下载路径：https://github.com/biocore/sortmerna/tree/master/data/rRNA_databases
3. 通过idba_tran v1.1.3 (https://github.com/loneknightpy/idba) 对转录组数据进行组装得到每个样本的转录组序列$Sample_assembly/contig.fa：
   # Installation
   $ conda install -c bioconda idba
   # Assembly
   $ $IDBA/bin/fq2fa –merge $Sample_rmRNA_R1.fq $Sample_rmRNA_R2.fq $Sample_merge.fa
   $ $IDBA/bin/idba_tran -r $Sample_merge.fa -o $Sample_assembly --pre_correction --mink 20 --maxk 60 --step 10 --num_threads 20
4. 组装得到的Contig中包含非mRNA信息，使用prodigal v2.6.3 (https://github.com/hyattpd/Prodigal) 对蛋白质编码区进行预测：
   # Installation
   $ conda install -c bioconda prodigal
   # ORFprediction
   $ prodigal -d $Sample_nul.fa -I $Sample_assembly/contig.fa -m -p meta
5. 预测到每个个体样品的基因序列后，利用CD-HIT v4.8.1 (http://weizhongli-lab.org/cd-hit/) 对所有样品的基因进行聚类构建非冗余基因集 (gene catalogue.fa)：
   # Installation
   $ conda install -c bioconda cd-hit
   # Gene clustering
   $ cat $Sample*_nul.fa > all_gene.fa
   $ cd-hit-est -I all_gene.fa -c 0.95 -aS 0.9 -M 0 -o all_gene_nr -T 40
   $ awk 'BEGIN{a=1}{if($0~/>/){print ">Unigene_"a;a+=1}else{print $0}}' all_gene_nr.fa >gene catalogue.fa
6. 样本中基因的表达量通过将reads比对到基因集 (gene catalogue.fa) 获得，通过将测序reads比对到基因序列上得到样品中基因表达量 (sample_fpkm.txt)。即使用RSEM v1.3.1 (https://github.com/deweylab/RSEM) 对bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) 的比对结果进行统计，获得每个样本中每个基因的 reads counts，同时计算FPKM。FPKM (全称expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) 是每百万 fragments 中来自某一基因每千碱基长度的fragments 数目，其同时考虑了测序深度和基因长度对fragments计数的影响，是对reads counts进行标准化处理的目前最为常用的基因表达水平估算方法 (Trapnell et al., 2010)。
   # Installation
   $ conda install -c bioconda rsem= 1.3.1
   # Expression calculate
   $ rsem-calculate-expression -p 70 --bowtie2 --paired-end \
   $Sample_rmRNA_R1.fq $Sample_rmRNA_R2.fq rsem.ref $Sample
   $ for i in `$Sample*.genes.results`; \
   do cut -f 1,7 $i >${i/gene.results/FPKM.txt}; \
   done
   $ python merge_metaphlan_tables.py $Sample*.FPKM.txt > all_sample_fpkm.txt
   (脚本下载链接: http://github.com/biobakery/MetaPhlAn/blob/master/metaphlan/merge_metaphlan_tables.py)
7. 使用diamond软件将非冗余的Unigenes 序列与 NCBI-NR 数据库进行比对 (设定阈值为e-value ≤ 1e-5)，采取最近公共祖先LCA (Lowest Common Ancestor)算法获得基因序列的物种注释信息：
   # Installation
   $ mkdir $HOME/gocode
   $ export GOPATH=$HOME/gocode/bin
   $ export PATH=$PATH:GOPATH
   $ go get github.com/emepyc/Blast2lca/blast2lca
   # KEGG annotation
   $ diamond blastp –threads 20 -q Unigenes_pro.fa -d nr_*version.dmnd -o Uni-genes_vs_nr_blt.txt --max-target-seqs 10 –evalue 1e-5 --outfmt 6 qseqid qlen sseqid stitle slen pident length mismatch gapopen qstart qend sstart send eval-ue bitscore
   $ blast2lca -input Unigenes_vs_nr_blt.txt -o Unigenes_vs_nr_blt.tax -ms 50 -me 0.000001 -g2t gi_taxid-March2015X.bin
   (gi号版本：https://pan.baidu.com/s/1Ekt8aSsoGjkZwr6PsZDKPA，提取码：i21t )
   (accession版本:https://pan.baidu.com/s/1WdZEaI1KPIMQ0PEv31lduw，提取码：5xfu)
8. 通过将基因序列和特定数据库 (如KEGG) 进行比对，完成基因功能注释。序列比对使用diamond v0.1.10 (https://github.com/bbuchfink/diamond) 软件进行。
   KEGG数据库 (https://pan.baidu.com/s/1jnulGNSQ3qDfoB3b76a0Dg，提取码：jzal,链接永久有效。
   # Installation
   $ conda install -c bioconda diamond=0.9.10
   # Annotation
   $ diamond makedb –in meta.pep -d meta.dmnd
   $ diamond blastp -d meta.dmnd -q Unigenes_pro.fa -k 1 -p 50 -f 6 -o Uni-genes_vs_kegg_blt.txt
   
   

三、结果分析
通过和KEGG 数据库中的KO数据库进行比对，得到不同层次的基因功能分类 (图2)。以一种酸性土壤样本为例，从KO level 1层次来看，土壤微生物的大部分活性基因与新陈代谢 (Metabolism) 相关，占总体基因表达量的45.7 %，其次为遗传信息加工 (Genetic information processing，8.5 %)，环境信息加工 (Environmental information processing，9.6 %) 和细胞过程 (Cellular processes，8.2 %)。KO level 2 层次上的基因功能分类如图2。 从图中可以看出，微生物新陈代谢相关的活动主要表现为能量代谢、碳水化合物代谢和氨基酸代谢；遗传信息加工主要表现在蛋白翻译；环境信息加工主要包括信号转导和膜转运；而细胞过程相关的基因主要参与了cellular community-eukaryotes和细胞迁移 (cell motility)。


图2.以酸性土样本为例的土壤微生物活性基因在KEGG level 2水平上的功能分类

四、土壤宏转录组研究的优点

1. 当与室内培养控制实验和高通量测序结合使用时，可以估算群落中特定微生物正在进行的积极的转录过程。
2. 排除了死亡微生物细胞残体，休眠体的影响。
3. 能够捕捉土壤特定类群内部动态变化。
4. 直接评估土壤微生物活性，包括对于干扰或者暴露等情况的响应。
   

注意事项

为防止DNA交叉污染，DNA污染物的去除很重要，纯化的RNA应该直接用PCR检测。 琼脂糖电泳缺乏检测复制片段表明缺乏检测到的交叉污染DNA。 如果检测到DNA，需要进一步使用DNase I分离RNA。

致谢

本实验得到国家自然科学基金项目（91951109）的资助。

参考文献

1. Mettel, C., Kim, Y., Shrestha, P.M. and Liesack, W. (2010). Extraction of mRNA from soil. Appl Environ Microbiol 76: 5995-6000.
2. Peano, C., Pietrelli, A., Consolandi, C., Rossi, E., Tagliabue, L., De Bellis, G.D. and Landini, P. (2013). An efficient rRNA removal method for RNA sequencing in GC­rich bacteria. Microb Inform Exp 3:1.
3. Trapnell, C., Williams, B.A., Pertea, G., Mortazavi, A.,Kwan, G., van Baren, M.J., Salzberg, S.L., Wold, B.J., and Pachter,L. (2010).Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol 28: 5.