1. 转录组数据质量控制分析的内容

主要包含几个方面：

1. 测序数据的质量控制，推荐使用 Trimmomatic 对低质量数据进行截除；
2. 是否是链特异性测序，以及链特异性测序的种类；
3. 转录组数据对参考基因组的匹配情况，包括分析匹配到基因内或基因间区的数据量；
4. 转录组数据中 rRNA 数据含量；
5. 转录子从 5' 到 3' 端的覆盖度分析；
6. 转录组测序数据量的饱和度分析。

本文主要讲解后面 4 个内容的分析方法。

2. 使用 RNA-SeQC

下载 RNA-SeQC 及其说明文档

$ wget http://www.broadinstitute.org/cancer/cga/tools/rnaseqc/RNA-SeQC_v1.1.8.jar
$ wget http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/RNA-SeQC_Help_v1.1.2.pdf

程序运行需求的输入文件：

genome.fasta        基因组的序列文件，同时需要有相应的 .fai 和 .dict 文件。
genome.gtf          基因组注释文件，该软件仅支持 GTF 格式。标准的 gtf 文件必须有 start_codon, stop_codon 和 CDS 内容，但本软件要求 gtf 文件必须还含有 exon 内容。
accepted_hits.bam   转录组数据比对到基因组的结果，例如 tophat 的结果。此 bam 文件必须含有 RG 参数。
rRNA.fasa           rRNA 序列，用于计算 rRNA 测序数据含量。此文件不是必须的。
rRNA_intervals.list rRNA 的结构信息文件，每行以 "chromosome_id:start-end" 信息表示 rRNA 结构。此文件不是必须的。但和上一个文件必须要二选一。

准备输入文件

得到基因组的 .fai 和 .dict 文件
$ samtools faidx genome.fasta
$ java -jar /opt/biosoft/picard-tools-1.95/CreateSequenceDictionary.jar R=genome.fasta O=genome.dict

若运行 tophat 没有输入 RG 参数，则运行下面命令得到 RG 参数
$ java -jar /opt/biosoft/picard-tools-1.95/AddOrReplaceReadGroups.jar I=accepted_hits.bam O=accepted_hits.RG.bam ID=A1 LB=A1 PL=ILLUMINA SM=A1 PU=run

若仅有 GFF3 文件，没有标准的 GTF 文件，需要自己编写程序进行转换
$ gff3ToGtf.pl genome.fasa genome.gff3 > genome.gtf
$ perl -p -i -e 's/^\s*$//' genome.gtf
若不去除 genome.gtf 文件中的空行，则会出现错误提示： No attributes were found in the GTF file on line 。

运行 RNA-SeQC 命令

对一个数据运行 RNA-SeQC，并结果输出到制定文件夹
$ java -jar RNA-SeQC_v1.1.8.jar -o A1 -r genome.fasta -s "A1|accepted_hits.RG.bam|A1_note" -t genome.gtf -BWArRNA rRNA.fasta

对多个数据运行 RNA-SeQC
$ echo -e "A1\tA1.bam\tA1_note
A2\tA2.bam\tA2_note
B1\tB1.bam\tB1_note
B2\tB2.bam\tB2_note" > samples.txt
$ java -jar RNA-SeQC_v1.1.8.jar -o rna-seqc_out -r genome.fasta -s samples.txt -t genome.gtf -BWArRNA rRNA.fasta

注意，多个数据运行 RNA-SeQC，输入的 samples.txt 文件内容为：
第一行信息为： "Sample ID\tBam File\tNotes", 这行头部会被忽略，仅起注释作用。
后面多行则为相关的信息，用3列表示，用tab分割，注意Bam文件必须为绝对路径。