【进化基因组学】比较转录组--数据预处理

前言

本文系视频的第四讲内容----数据前处理。

在开始之前，大家可以想一想如果是在sanger测序的时候，我们第一步就是需要确定样品，然后是对样品提取DNA，检验DNA质量，之后便是根据引物去配置pcr体系，然后跑pcr，跑完之后，我们还会进行跑胶检验，看看是否真的p出条带了。假如一切顺利，p出了条带，那我们就会拿去测序公司对样品进行测序。测完序之后，我们会放进读峰图的软件，然后把测序质量可靠的保留下来，作为后续分析的基础。

其实我们做数据分析也是一样的思维：首先是获得一份初始数据（即我们提取了样品的DNA），经过数据初步质检（DNA质量），然后经过数据清洗（我们设计了引物，只p出我们想要的序列片段），再检验一下数据质量（跑胶），确定初始数据可用之后（存在目标条带），会经过一个上游分析（读峰图），得到我们下游分析所需要的数据（切除测序质量较差的碱基）。

而对应到我们的比较转录组，自然就是先经过测序得到原始数据，又或者是使用已经发表的转录组数据（初始数据），经过数据初步质检（fastqc确定碱基质量），切除接头与低质量碱基（数据清洗），然后再检验一次质量（fastqc确定质量是否满足后续分析），经过序列拼接得到转录本，而后再去冗余，再经过编码区预测（上游分析），得到后续进化分析所需要的数据---转录本（unigene），编码序列（cds），蛋白序列（pep）

原始数据获取

若是正常的课题测新物种，一般都是交由测序公司进行RNA提取与建库测序。至于样品的要求，做比较转录组对时期不是那么严格，只要能够获取大量生物体RNA即可。

这样说可能比较懵逼，就拿个例子，假如我现在是出差去野外采样，我采了一两株植物，然后他们带回来实验室已经差不多被折腾得半死了，但依旧还是活着的，这样我就还是可以直接拿去提取RNA，进行比较转录组的分析。

而这一次我们的测试数据选择来自NCBI的两个物种Arabidopsis halleri与Arabidopsis lyrata，这两个物种其实都有基因组公布了，各位如果有兴趣可以拿比较转录组的结果跟已发表基因组的结果做对比，就可以知道比较转录组的精确度如何了。

# SRR12062107(Arabidopsis halleri)
https://sra-download.ncbi.nlm.nih.gov/traces/sra9/SRR/011779/SRR12062107
# SRR12062120(Arabidopsis lyrata)
https://sra-download.ncbi.nlm.nih.gov/traces/sra74/SRR/011779/SRR12062120

下载工具各位可以任选，反正我用迅雷就挺好的。

SRA转fastq

我们从ncbi上下载得到的两个测序数据还不能直接用，因为它们并不是我们常用的测序数据格式，而是另一种格式SRA---ncbi的通用数据存储格式，采用二进制压缩，而我们的fastq文件则是文本类的文件，故而我们需要借助一些软件对数据类型进行转换。

NCBI官方提供的一种软件叫做sratoolkit，它里面内置的fastq-dump模块就可以将sra文件转换成为fastq文件，具体下载地址如下：

https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.8.2/sratoolkit.2.8.2-centos_linux64.tar.gz

安装也是非常简单，直接解压就可以用了，而我们为了方便调用，通常需要添加进全路径

#解压
tar xzvfsratoolkit.2.8.2-centos_linux64.tar.gz
#添加进全局变量
echo "exportPATH=$PATH:/Biosoft/sratoolkit.2.8.2-centos_linux64/bin" >>~/.bash_profile
source ~/.bash_profile

但这个软件写的足够落后，单线程，效率还低，一个sra文件需要“山无棱，天地合”的时间才刚给你转成fastq文件，故而我们需要给它打一管鸡血--pfastq-dump

# 下载地址
https://github.com/inutano/pfastq-dump
# 下载完之后，将bin文件夹的pfastq-dump放入sratoolkit的bin文件夹，并修改执行权限即可
chmod a+x pfastq-dump

进行数据转换

pfastq-dump  SRR12062107 --split-3 -t 4 -O ./
# 如果是单端测序则不需要添加--split-3
# -t表示线程数，本次使用了4个线程
# -O输出目录，./表示当前目录

初步质检

首先，我们先来看一份质量报告

一开始你可能会觉得很神奇，怎么这个质检软件可以得到这么多信息，那么接下来我们就来揭一下质检软件神秘的面纱。 有因必有果，你的报应就是我，要了解质检软件是怎么进行统计的，我们就需要了解一下fastq格式是怎么一回事

ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG
+
FFKKKFKKFKF<KK<F,AFKKKKK7FFK77<FKK,<F7K,,7AF<FF7FKK7AA,7<FA,,
# 第一行就是测序的坐标信息，即告诉你这条reads的名字是什么
# 第二列就是我们测到的序列
# 第三列就一个加号，附加信息，正常没有
# 第四列，质量信息，对应着上面各个碱基

质检软件其实就是通过统计这几个信息进行，每条序列的长度是固定的，每个碱基的质量也是固定的，其实软件只要做几个累加数值就可以得到一份质检报告。

那么有同学还有问题，那这个接头序列是怎么看的，其实就是高通量测序的时候，我们会加入接头序列，而这些序列通常会内置到质检软件里面，软件只要发现了一样的序列在两端，即可进行统计，自然也就知道数量了。

既然已经知道原理了，那么我们要怎么样去得到这样一份质检报告呢？很简单，就是我们的老牌质检软件---fastqc。

# 安装，使用conda安装即可
conda install -c bioconda fastqc
# 使用
fastqc -o ./ -t 6 SRR12062107_1.fastq.gz  SRR12062107_2.fastq.gz
# -o质检报告输出目录
# -t 指定线程数

数据清洗

由质检报告我们可以知道，数据中存在相当部分低质量碱基，还有很少的接头序列，为了不影响后续分析的准确度，我们需要把这些“拖后腿”的序列全部剔除 可以进行序列过滤的软件很多，大部分都可以从conda进行安装，本次进行示范的是trimmomatic，一个java软件，下载即用，无需编译

# 下载地址：
http://www.usadellab.org/cms/index.php?page=trimmomatic
            

解压之后是一个jar文件跟一个adapters的文件夹（内置了常用接头序列），接着我们来了解一下trimmomatic的常用参数

# ILLUMINACLIP: 过滤 reads 中的Illumina 测序接头和引物序列，并决定是否去除反向互补的 R1/R2 中的 R2
# SLIDINGWINDOW: 从 reads 的 5'端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗
# MINLEN: 如果经过剪切后 reads的长度低于阈值则丢弃这条 reads
# AVGQUAL: 如果 reads的平均碱基质量值低于阈值则丢弃这条 reads

这里关于引物的选择，我们可以来看看trimmomatic的内置文件夹 PE是双端测序，SE是单端测序，此处应该不难知道意思，而最主要是三种类型的接头要怎么选，其实就看看fastqc那一项adapter的报告即可，正常illumina的在TruSeq里面任选就行。

下面看看本次的第一个批量执行的脚本

for fn in *_1.fastq.gz # 读入所有1端数据
do
  sample=${fn%%_*} # 将全传入参数截取_前部分传入，即SRR号
  java -jar/public1/users/houzw/Trimmomatic/Trimmomatic-0.38/trimmomatic-0.38.jar \
       PE \
       -threads 4 \
       -basein \
         ${sample}_1.fastq.gz \
         ${sample}_2.fastq.gz \
       -baseout \
        ${sample}.R1.clean.fq.gz\
        ${sample}.R1.fail.fq.gz\
        ${sample}.R2.clean.fq.gz \
        ${sample}.R2.fail.fq.gz\
      ILLUMINACLIP:/public1/users/houzw/Trimmomatic/Trimmomatic-0.38/adapters/NexteraPE-PE.fa:2:30:10\ # 全路径指定引物与切除步长
       MINLEN:50 \
       AVGQUAL:20
done

等到清洗完成再投入fastqc进行一次质检，满足要求了即可进入下一步分析

序列拼接（本节配图有一处错误，希望仔细看的同学可以看出来）

怎么把序列拼接做个类比，个人觉得是拼乐高玩具比较合适。

有参拼接相当于给了你一份图纸，让你知道每一步应该怎么走，而无参拼接则像是把图纸弄丢了，然后要拼出来，怎么拼纯靠逻辑推理，中间可能拼的顺序是看似正确的，还有些部件就完全联系不上了。

拼接的原理其实说得粗暴一点就是：谁像跟谁连。以无参转录组的老牌拼接软件：Trinity为例，稍微详细一点就是下面这样：

•Inchworm:从reads中提取所有的重叠k-mers，根据丰度递减的顺序检查每个k-mers，然后将重叠的k-mers延长到不能再延长，称为一个contig

•Chrysalis: 将上一部生成的contig聚类，对每个类构建deBrujin graph

•Butterfly: 根据构建的deBrujin graph ，寻找具有可变剪接的全长转录本，同时将旁系基因的转录本分开

原理说完，我们就来看看trinity正式使用应该怎么做吧

# trinity的安装有两种方式，一种是下载源码进行编译（不太推荐）
https://github.com/trinityrnaseq/trinityrnaseq

# 第二种是使用conda创建虚拟环境，
# 主要是由于trinity版本迭代，
# 对于perl的版本要求可能跟主环境不同，故而需要单独创建一个环境，避免环境交叉污染
conda create -n trinity
conda activate trinity
conda install -c bioconda trinity

# trinity的常用参数主要是这几个
--seqType 指定类型是fa还是fq
--max_memory 指定内存，以G为单位
--left R1端输入文件
--right R2端输入文件
--SS_lib_type 链特异性建库时选择的参数
--CPU 指定线程
--output 指定输出目录

各位如果在一些老旧的服务器上运行trinity，切记要跑一个多线程跟一个单线程的，多半由于I/O性能不行，多线程会空转，所以推荐跑两种模式，稳一稳没什么不好

for fn in *.R1.clean.fq.gz
do
  Trinity --seqType fq \
        --left $fn \
        --right ${fn%%.*}.R2.clean.fq.gz\
        --CPU 4 \
        --max_memory 10G
done

转录本去冗余

首先，我再次强调这个思维：转录本其实就是相当于一个没有了间隔区的“简化”基因组。

但我们拼接得到的trinity结果却还包含着一个基因的许多个可变剪切，我们的目的并不是研究可变剪切，而是研究基因水平的序列变化，故而多个可变剪切也只能算是一个基因。

我们一般认为可变剪切里最长的序列是最完整反应整个基因的序列，故而去冗余的思想其实就是保留最长的可变剪切转录本。

trinity内置就有去重脚本，我们也可以根据trinity的命名原则去去重，也可以通过第三方的软件进行去重。

Trinity的命名规则
TRINITY_DN1000_c115_g5_i1
TRINITY_DN1000_C115_g5是一个基因
i1，i2....表示该基因的不同转录本
思路：保留最长转录本，即在同一基因的不同转录本找最长的即可，
可以用字典key先保存基因名，value为序列；
每次比较的时候，用len()比较不同转录本的大小，
大则留下，小则剔除，遍历完整个文件就可以依次输出

软件去重的话，目前用得最多的是cd-hit

# 二进制版，解压即用
https://github.com/weizhongli/cdhit/releases/download/V4.6.7/cd-hit-v4.6.7-2017-0501-Linux-binary.tar.gz
# 基本命令参数：
-i输入文件，必须为fasta格式
-o输出文件路径
-c相似度
-n两两序列比对时选择的word size
-d 0时表示使用fasta标题中第一个空格前的字段作为序列名字
-M 内存，以m为单位
-T 线程数

从我个人的使用情况来说，一般无参转录本拼接出来，相似度选择过滤在85%-90%是一个比较合适的阈值，而wordsize 一般选5-7都可以

cd-hit -i input.fa -o input.db85.fa -c 0.85 -n 5 -M 4000 -d 0 -T 4

拼装结果评估

对于很多人来说，总是记得对原始数据进行质控，却在后面的环节把这一步忘得一干二净，其实质控是一个重要的环节，也是一个必须有的习惯，这样才能对数据的变化情况有所把握，而不是简简单单的在跑流程。

对于无参转录组的拼装评估，其实最简单的就是拿去完冗余的转录本---Unigene作为reference transcript，然后使用clean data进行回帖，统计mapping rate，然后再统计N50，N90等指标即可，通常mapping rate可以达到90%即为较优秀的拼接结果

预测ORF

在得到一个转录本之后，我们就相当于得到了一整个物种的基因信息，而物种进行生命活动需要进行翻译，将遗传信息转化为表达载体，即从核DNA--->cds--->pep的过程。

现在我们已经相当于拥有核DNA的信息，需要我们按照表达规则，将cds与pep全部翻译出来，这个时候我们会用到一个叫做transdecoder的软件对整个转录本进行预测，通常这一步会结合注释结果再进行整合。

# 安装
conda install -c bioconda transdecoder
# 使用LongOrf模块进行预测
TransDecoder.LongOrfs -t unigene.fa
# 运行结果会生成两个文件夹
# 在unigene.fa.transdecoder_dir里面存放着预测的cds,pep,gff文件

总结

本节课，我们通过一个简单的转录组到转录本的流程，对整个比较转录组预备分析阶段做了一个初步的了解，贯穿始终的灵魂就是质量检测，这个是我们每一步都不能掉以轻心的，只有牢牢把握数据的质量，我们才能更好的了解和完成整个项目，而不是一个代码流程工。