做驴转录组数据然后脑袋被驴踢了搞出来几万个差异

下面实习生投稿

这几天在复现一篇文章《Single-Cell RNA-Seq Revealed the Gene Expression Pattern during the In Vitro Maturation of Donkey Oocytes》，在对数据完成了过滤、比对和定量后，开始进行下游分析。

在我样本检测一顿输出后，拿到以下五个图，感觉数据还不错

问题出现

但当我开始做差异分析的时候，问题就出现了，在我将FoldChange的阈值设置为2，pvalue设置为0.01时，上调的基因有7千多个，下调的基因也有六千五百多个，尽管上下调基因数量和normal基因数量比例还算合适，但这个数量也太离谱了。

我开始往前寻找问题的根源，在刚读取表达矩阵时，我的矩阵居然有57万行，但是对于当时的我来说，尚未发现异常。

因为本来就是会随后对表达矩阵进行了预处理，过滤掉了至少在75%的样本中都不表达的基因

keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
filter_count <- rawcount[keep,]
filter_count[1:4,1:4]

此时再查看，剩下4万行左右，但这个数量还是超纲了。（正常情况下应该是2万个基因，不过主要是取决于gtf文件的记录情况）

查看了一下表达矩阵，嘶，这些居然是外显子......

完蛋，做成了差异的外显子了。（非常懊恼啊，简直是脑子被驴踢了！）

罪源

究其原因，还是在上游分析定量时的这条命令 ：

nohup featureCounts -T 20 -g 'ID' -p -a /home/xiaowang/donkey_oocyte/1_genome/genomic.gff -o counts.txt /home/xiaowang/donkey_oocyte/4_mapping/*.sam >count.log 2>&1 &

-g 'ID'这个参数直接导致结果的行名为exon，而不是基因名。换一句话说，featureCounts默认定量的水平为Meta-features。

feature指的是基因组区间的最小单位，比如exon; 而meta-feature可以看做是许多的feature构成的区间，比如属于同一个gene的外显子的组合，也可以是不同转录本。

• -t "exon" -g "ID" 相当于手动将定量的level又降到exon了
• -t "exon" -g "gene_id" 则是将定量到的exon汇总到gene水平上

解决办法

既然如此，提供两条解决办法

1. 将相同基因的外显子的counts合并
2. 重新定量 此时我的命令为

nohup featureCounts -T 20 -t "exon" -g 'gene_id' -p -a /home/xiaowang/donkey_oocyte/1_genome/genomic.gtf -o counts.txt /home/xiaowang/donkey_oocyte/4_mapping/*.sam >count.log 2>&1 &

过滤后的表达矩阵

上下调基因数量也就因此降了下来。

fc_cutoff<-2 pvalue<-0.01

以上。

另外，值得一提的是实习生他也有自己的公众号，有一个保研专栏，感兴趣可以去看看！