提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  数据库

当前栏目

测序数据GC含量异常该如何处理？

异常 数据
2023-03-07 09:16:21 时间

我们在对测序数据进行fastqc质控时，会比较关注样本的GC含量，较好的数据如下图所示

Normal data

事实上，我们的测序结果会受到建库方式、样本质量等很多因素的影响，其GC含量也会千差万别。下图是真实测序结果的fastqc报告，很明显其GC含量有多个峰值，相当的“红”。我们在进行后续分析之前，首先要对其进行一定处理，去掉杂峰的影响。

Abnormal data

为了便于测试，我选取了一个样本的双端测序结果进行演示，如果想自己尝试可翻到文末数据获取部分下载数据。

rawdata_qc

下图是双端测序结果的fastqc质控报告中的GC含量部分，我们可以看到有几个异常峰值。

曾老师曾经在一篇文章中提到右峰有可能是rRNA，就这个思路我们不妨将测序数据比对到不同的RNA上，没有比对上的数据即为去掉了对应RNA的数据。

需要使用到hisat2的--un-conc-gz参数~

RNA数据的下载

由于本次我们使用的是人的数据，所以在Nucleotide中输入Homo sapiens，如果是其他物种如拟南芥输入Arabidopsis thaliana（对应学名）即可。

将需要的RNA选中类型选中

以rRNA为例，输出为FASTA格式，再传到服务器按照常规比对方法比对即可。

rm_rRNA_qc

首先我们比对到rRNA上，比对的代码和结果如下，共有14.69%的reads比对到了rRNA上

# 建立索引
hisat2-build -p 16 rRNA.fasta rRNA

# 比对
hisat2 -x /home/xiaowang/GSE145894/4.rRNA_rm/genome_index/rRNA -1 /home/xiaowang/GSE145894/2.raw_fq/SRR11178348_1.fastq.gz -2 /home/xiaowang/GSE145894/2.raw_fq/SRR11178348_2.fastq.gz --un-conc-gz ./test_rm.fq.gz -p 16 -S ./test.sam

# 比对结果
49273485 reads; of these:
  49273485 (100.00%) were paired; of these:
    42251347 (85.75%) aligned concordantly 0 times
    632294 (1.28%) aligned concordantly exactly 1 time
    6389844 (12.97%) aligned concordantly >1 times
    ----
    42251347 pairs aligned concordantly 0 times; of these:
      333 (0.00%) aligned discordantly 1 time
    ----
    42251014 pairs aligned 0 times concordantly or discordantly; of these:
      84502028 mates make up the pairs; of these:
        84072885 (99.49%) aligned 0 times
        33143 (0.04%) aligned exactly 1 time
        396000 (0.47%) aligned >1 times
14.69% overall alignment rate

# 质控
fastqc test_rm.fq*

以下是同一样本双端数据的fastqc结果

通过与原始数据对比，我们已经可以发现右峰被去除，说明右峰为rRNA。

rm_rRNA_ncRNA

在此基础上，我们再将数据比对到ncRNA上，共有11.19%reads比对到ncRNA上。

# 比对结果
42251347 reads; of these:
  42251347 (100.00%) were paired; of these:
    37862094 (89.61%) aligned concordantly 0 times
    1965637 (4.65%) aligned concordantly exactly 1 time
    2423616 (5.74%) aligned concordantly >1 times
    ----
    37862094 pairs aligned concordantly 0 times; of these:
      24274 (0.06%) aligned discordantly 1 time
    ----
    37837820 pairs aligned 0 times concordantly or discordantly; of these:
      75675640 mates make up the pairs; of these:
        75050031 (99.17%) aligned 0 times
        400055 (0.53%) aligned exactly 1 time
        225554 (0.30%) aligned >1 times
11.19% overall alignment rate

此时去掉了rRNA和ncRNA的fastqc结果已经符合要求了。

rm_ncRNA

以下为原始数据比对到ncRNA后的fastqc结果，中间的小峰被去除。

同时我又尝试将测序数据直接比对到tRNA和cRNA上，比对率都为0，说明序列中都不包含tRNA和cRNA。

需要注意的是，生物体中的RNA理论上按照是否编码蛋白质将其分为编码RNA（coding RNA）和非编码RNA（non-coding RNA, ncRNA）两大类，前者指mRNA，后者则包括很多种类，如众所周知的tRNA和rRNA，参与mRNA剪接的snRNA，参与RNA修饰的snoRNA等。非编码RNA种类很多，按照长度可分为两大类：大于200 nt（核苷酸）的称为长非编码RNA（lncRNA），小于200 nt的称为小非编码RNA（small ncRNA），50 nt以下的还可称为tiny ncRNA，如siRNA、miRNA和piRNA等。

但我在测试的时候发现，上文中提到的rRNA、ncRNA、tRNA、cRNA在数据上没有交集，也就是说将去掉ncRNA后的数据比对到tRNA和cRNA的比对率为0并不是因为ncRNA的数据包含tRNA和cRNA，而是测序数据本身。

数据获取

以上演示数据来自GEO数据集GSE145894，按照以下方法即可下载

拿到Accession List以后进入服务器通过prefetch下载

cat SRR_Acc_List.txt | while read id
do
 echo prefetch ${id} -O ./
done

参考资料：

  1. RNA-seq数据分析完全指北-03：去除奇怪的RNA
  2. RNA-seq的fastq文件里面为什么有gc含量的双峰

相关文章