前言
本文系视频的第四讲内容----数据前处理。
在开始之前，大家可以想一想如果是在sanger测序的时候，我们第一步就是需要确定样品，然后是对样品提取DNA，检验DNA质量，之后便是根据引物去配置pcr体系，然后跑pcr，跑完之后，我们还会进行跑胶检验，看看是否真的p出条带了。假如一切顺利，p出了条带，那我们就会拿去测序公司对样品进行测序。测完序之后，我们会放进读峰图的软件，然后把测序质量可靠的保留下来，作为后续分析的基础。
其实我们做数据分析也是一样的思维：首先是获得一份初始数据（即我们提取了样品的DNA），经过数据初步质检（DNA质量），然后经过数据清洗（我们设计了引物，只p出我们想要的序列片段），再检验一下数据质量（跑胶），确定初始数据可用之后（存在目标条带），会经过一个上游分析（读峰图），得到我们下游分析所需要的数据（切除测序质量较差的碱基）。
而对应到我们的比较转录组，自然就是先经过测序得到原始数据，又或者是使用已经发表的转录组数据（初始数据），经过数据初步质检（fastqc确定碱基质量），切除接头与低质量碱基（数据清洗），然后再检验一次质量（fastqc确定质量是否满足后续分析），经过序列拼接得到转录本，而后再去冗余，再经过编码区预测（上游分析），得到后续进化分析所需要的数据---转录本（unigene），编码序列（cds），蛋白序列（pep）
原始数据获取
# 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
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
# 下载地址
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输出目录，./表示当前目录
初步质检
首先，我们先来看一份质量报告
ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG
+
FFKKKFKKFKF<KK<F,AFKKKKK7FFK77<FKK,<F7K,,7AF<FF7FKK7AA,7<FA,,
# 第一行就是测序的坐标信息，即告诉你这条reads的名字是什么
# 第二列就是我们测到的序列
# 第三列就一个加号，附加信息，正常没有
# 第四列，质量信息，对应着上面各个碱基
# 安装，使用conda安装即可
conda install -c bioconda fastqc
# 使用
fastqc -o ./ -t 6 SRR12062107_1.fastq.gz SRR12062107_2.fastq.gz
# -o质检报告输出目录
# -t 指定线程数
数据清洗
# 下载地址：
http://www.usadellab.org/cms/index.php?page=trimmomatic
# ILLUMINACLIP: 过滤 reads 中的Illumina 测序接头和引物序列，并决定是否去除反向互补的 R1/R2 中的 R2
# SLIDINGWINDOW: 从 reads 的 5'端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗
# MINLEN: 如果经过剪切后 reads的长度低于阈值则丢弃这条 reads
# AVGQUAL: 如果 reads的平均碱基质量值低于阈值则丢弃这条 reads
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_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
拼装结果评估
预测ORF
# 安装
conda install -c bioconda transdecoder
# 使用LongOrf模块进行预测
TransDecoder.LongOrfs -t unigene.fa
# 运行结果会生成两个文件夹
# 在unigene.fa.transdecoder_dir里面存放着预测的cds,pep,gff文件
总结
本节课，我们通过一个简单的转录组到转录本的流程，对整个比较转录组预备分析阶段做了一个初步的了解，贯穿始终的灵魂就是质量检测，这个是我们每一步都不能掉以轻心的，只有牢牢把握数据的质量，我们才能更好的了解和完成整个项目，而不是一个代码流程工。