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文件必须为绝对路径。