无参转录组之转录本质控
刘小泽写于19.2.2
今天主要看看转录本拼接质量评估
前言
我们利用trinity简单的命令拼接好转录本后，想要知道我们拼接的质量如何，如果不好可能需要寻找原因，重新调整参数再重新拼接。这一步是很关键的，因为拼接的转录本相当于有参里面的参考转录组，试想如果参考都做不好，那岂不是上梁不正下梁歪？
首先我们可以看看组装了多少条转录本
grep '>' $wkd/assembly/trinity_out_dir/Trinity.fasta | wc -l
但这仅仅是最粗略的方式，因为随着测序深度的加大，得到的contigs越多，因此可以拼接更多的转录本
另外可以看看基本的统计值
TrinityStats.pl $wkd/assembly/trinity_out_dir/Trinity.fasta
第一部分Counts of transcripts, etc. 结果中的 Total trinity transcripts 和上面那个命令结果一致【正常数据中转录本数量小于20w是正常的，如果数量达到了30w、40w条，就需要先用corset软件进行聚类】
结果中还有组装的gene数 Total trinity 'genes' ，可以看到transcripts的数量比genes的数量多，因为存在一个基因的可变剪切（这种情况在昆虫和哺乳动物中比较常见）
另外还有第二部分内容Stats based on ALL transcript contigs: N50值（at least half of the assembled bases are in contigs of at least that contig length
累加后长度超过转录组总长度一半的contig的长度就是N50）
【正常情况下，N50应该是1k左右】
说到评估，我们一般会想到：和其他拼接工具的结果进行对比（就好像在有参中使用3-5中工具进行比对，看看分别的比对率），或者使用不同的参数看看结果异同（前提是对参数的设置比较了解）
总的来说，上面的方法得到的结果不是特别有用，有下面几种方法可以更有效地评价转录本质量
第一种 align
将reads重新比对到转录本上，看看比对率
一般来讲，至少应该有80%的原始数据可以在拼接的转录本中找到（经验值），剩下的没有拼接上的序列可能由于低表达导致没有足够的覆盖度进行拼接或者序列质量较低或者重复reads。
比对过程中，除了正确比对PE reads，还有可能仅仅一条read比对上，或者由于序列太短比对不上，如果利用bowtie、STAR进行比对，仅仅会计算成功比对上的PE reads。这里为了客观评价，我们需要知道所有的比对情况（包括单端比对和未比对上），可以利用Bowtie2
首先为转录本构建索引
bowtie2-build Trinity.fasta Trinity.fasta
然后进行双端比对，输出的sam转为bam
bowtie2 -p 10 -q --no-unal -k 20 -x Trinity.fasta -1 reads_1.fq -2 reads_2.fq \
2>align_stats.txt| samtools view -@10 -Sb -o bowtie2.bam
#-p 线程数
#-q input文件为fq(默认)
#--no-unal输出的sam文件中不记录未必对的reads信息
#-k 指定每条read匹配上的前20个记录，按匹配分值降序排列（descending order by alignment score）this mode can be effective and fast when user cares more about whether a read aligns (or aligns a certain number of times) than where exactly it originated.
bowtie2说明书：http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml
可以将比对结果放入IGV中查看，但首先需要排序、构建索引
ls $wkd/evaluate/1-bowtie2/*.bam | while read i; do
file=`basename $i`
#echo $file
surname=${file%%.bowtie2.bam}
#echo $surname
samtools sort $i -o $wkd/evaluate/1-bowtie2/${surname}.sorted.bam
samtools index $wkd/evaluate/1-bowtie2/${surname}.sorted.bam
done
samtools faidx $wkd/assembly/trinity_out_dir/Trinity.fasta
IGV
上面因为是选取的10000条序列，因此比对的结果一般，实际的数据更应该接近于下图
实际IGV
第二种 uniprot
将拼接转录本与已知蛋白序列库进行比对
核心思想就是检查拼接好的全长或者接近全长的转录本序列数目（因为这样的序列才有比对结果）。这种方法如果用在有参考转录本的模式物种（如人、小鼠）是很简单的，直接比对到参考转录组然后检查比对上的序列长度区间；非模式物种可以比对到近源物种的参考转录组上进行查看；更为一般的做法，就是将拼接的转录本比对到已知的蛋白数据库，看看能比对到多少蛋白序列，从而推断拼接质量
利用blast+和蛋白数据库SwissProt、TrEMBL 【关于这两个数据库：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC102476/】
简而言之，
SWISS-PROT: curated, high level of annotation, minimal level of redundancy ;
TrEMBL: computer-annotated upplement to SWISS-PROT, derived from the translation of all coding sequences (CDSs) in the EMBL except the CDSs already included in SWISS-PROT
#Build a blastable database:
makeblastdb -in uniprot_sprot.fasta -dbtype prot
# Perform the blast search, reporting only the top alignment:
cp $wkd/assembly/trinity_out_dir/Trinity.fasta ./
blastx -query Trinity.fasta -db uniprot_sprot.fasta -out blastx.outfmt6 \
-evalue 1e-20 -num_threads 6 -outfmt 6
# -max_target_seqs 1 这个参数可以加上试试 Maximum number of aligned sequences to keep
关于blast参数(- max_target_seqs)的设置还专门有一篇文章讲重要性：https://academic.oup.com/bioinformatics/advance-article-abstract/doi/10.1093/bioinformatics/bty833/5106166
结论就是：这个max_target_seqs参数有bug，不懂就用默认
检查结果
distribution of percent length coverage for the top matching database entries就是比对蛋白结果的序列长度分布，这么说可能有点拗口，简单理解就是：最左侧一列是蛋白序列比对到Trinity拼接转录本的覆盖度，第2列是落在某个覆盖度范围内的蛋白序列个数，第三列是一个蛋白数目累加值
analyze_blastPlus_topHit_coverage.pl blastx.outfmt6 Trinity.fasta uniprot_sprot.fasta
# 参数设置： blast+.outfmt6.txt query.fasta search_db.fasta
cat blastx.outfmt6.hist
#hit_pct_cov_bin count_in_bin >bin_below
100 71 71
90 22 93
80 11 104
70 11 115
60 19 134
50 22 156
40 38 194
30 25 219
20 62 281
10 14 295
如果比对的目标蛋白序列同时匹配多个Trinity转录本序列，并且都是best hits，那么结果也只计作一次作为highest BLAST bit score和longest match length
（Only the single best matching Trinity transcript is reported for each top matching database entry）
上面的结果主要理解成：
- 有22个蛋白，它们的序列的80%~90%可以比对到Trinity拼接的转录本
- 有93个蛋白，它们有超过80%的序列覆盖度（即可以比对到Trinity拼接的转录本），可以认为有93个比对到了全长转录本
- 有71个蛋白，它们超过90%的序列覆盖度
同样，可以进行Nr, Pfam，Nt数据库的比对，只要下载数据库=》建库=》比对，当然，如果下载的数据库（如：NCBI的nt、nr）已经构建好了索引，就可以跳过建库的步骤
统计升级
blast的过程存在这种情况：一条转录本比对到一个蛋白序列，另外还有几个不连续的片段（即：multiple high scoring segment pairs (HSPs)）。上面的blast流程只考虑了single best alignment这种情况，我们可以将每个转录本的HSPs组合在一起，然后基于此进行覆盖度计算
# Group the multiple HSPs
$TRINITY_HOME/util/misc/blast_outfmt6_group_segments.pl blastx.outfmt6 \
Trinity.fasta uniprot_sprot.fasta > blastx.outfmt6.grouped
# compute the percent coverage by length histogram
blast_outfmt6_group_segments.tophit_coverage.pl blastx.outfmt6.grouped
#hit_pct_cov_bin count_in_bin >bin_below
100 72 72
90 22 94
80 11 105
70 11 116
60 19 135
50 23 158
40 37 195
30 24 219
20 61 280
10 14 294
虽然得到的结果和上面单独blast差不多，但是目前是基于grouped HSPs/pair而非single HSP/pair
第三种 BUSCO
在近缘物种中，总有一些基因是保守的，BUSCO 软件根据OrthoDB 数据库，构建了几个大的进化分支的单拷贝基因集。将转录本拼接结果与该基因集进行比较，根据比对上的比例、完整性，来评价拼接结果的准确性和完整性。
安装：conda install busco
目前支持的数据库列表如下，选择研究物种的数据库进行下载。不同的数据库，评价结果可能完全不同。
BUSCO
mkdir $wkd/evaluate/5-busco && cd $wkd/evaluate/5-busco
# 下载数据库
wget http://busco.ezlab.org/v2/datasets/fungi_odb9.tar.gz
# 提取每个基因的最长转录本
get_longest_isoform_seq_per_trinity_gene.pl $wkd/assembly/trinity_out_dir/Trinity.fasta \
>longest_isoform.fasta
# 运行 BUSCO
# -l/--lineage: 数据库的路径
# -m:运行模式，geno|tran|prot
busco -i longest_isoform.fasta \
-l euarchontoglires_odb9 \
-o busco \
-m tran \
--cpu 2
结果主要看统计文件run_busco/short_summary_busco.txt，这里尽量选用真实数据进行测试，数据量太少得到的结果不准确
结果类似于：
# Summarized benchmarking in BUSCO notation for file longest_isoform.fasta
# BUSCO was run in mode: tran
C:90.0%[S:47.9%,D:38.2%],F:3.4%,M:6.7%,n:439
383 Complete BUSCOs (C)
216 Complete and single-copy BUSCOs (S)
167 Complete and duplicated BUSCOs (D)
20 Fragmented BUSCOs (F)
27 Missing BUSCOs (M)
430 Total BUSCO groups searched
C值(Complete BUSCOs)表示和BUSCO集相比的完整度，一般能达到80%以上。但是如果达不到也不需要担心，可以再多组装几个版本，选个最优解
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com
Welcome to our bioinfoplanet!
