对转录组数据进行质量控制分析

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

对转录组数据进行质量控制分析》上有1条评论

发表评论

电子邮件地址不会被公开。 必填项已用*标注