猪的单细胞分析如何过滤线粒体基因

当我们进行非模式物种的单细胞分析时，难免会遇到各种问题，其中，基因symbol命名不规范也将会给我们后续注释带来不小的麻烦。

比如我们在做猪的单细胞数据分析时，首先进行质控，当我们想要过滤线粒体基因时，采用常规的MT开头的基因进行过滤，会发现找不到基因。举个例子：

mito_genes=rownames(sce.all)[grep("^MT-", rownames(sce.all))] 
mito_genes #13个线粒体基因
# character(0)

其实，线粒体基因是真实存在的，只不过采用的命名方式不是MT开头，所以上面的代码失效了。

下面我们来演示一下在猪的单细胞分析如何过滤线粒体基因。

首先我们在猪的Ensembl数据库

数据库截图

发现存在一个MT线粒体gff3文件，希望在里面能找到些信息，解压并打开文件，果不其然，文件中的gene命名是以ENSSSCG开头的；

找一个基因测试下ENSSSCG00000018060

> mito_genes=rownames(sce.all)[grep("ENSSSCG00000018060", rownames(sce.all))] 
> mito_genes
# [1] "ENSSSCG00000018060"

真的存在这个基因，接下来就是如何从gff3文件中提取geneID，运行下列代码提取geneID：

grep -v "#" Sus_scrofa.Sscrofa11.1.105.chromosome.MT.gff3 |awk -F "[\t=:;]" 'BEGIN{OFS="\t"}$3=="gene"{print $1,$4,$5,$11}' >PIG.MT.gene.nemes.txt

MT RefSeq 3922 4876 ENSSSCG00000018065 MT RefSeq 5087 6128 ENSSSCG00000018069 MT RefSeq 6511 8055 ENSSSCG00000018075 MT RefSeq 8203 8890 ENSSSCG00000018078 MT RefSeq 8959 9162 ENSSSCG00000018080 MT RefSeq 9120 9800 ENSSSCG00000018081 MT RefSeq 9800 10583 ENSSSCG00000018082 MT RefSeq 10653 10998 ENSSSCG00000018084 MT RefSeq 11069 11365 ENSSSCG00000018086 MT RefSeq 11359 12736 ENSSSCG00000018087 MT RefSeq 12935 14755 ENSSSCG00000018091 MT RefSeq 14739 15266 ENSSSCG00000018092 MT RefSeq 15342 16481 ENSSSCG00000018094

所以，我们将用着13个线粒体基因ID。这个代码来自于我们《生信技能树》的单细胞数据分析工程师 林枫，欢迎大家前来咨询我们团队的单细胞数据分析服务哦

接下来我们进行实战

前面的PercentageFeatureSet()函数有两种指定模式：

• (1)pattern通过正则表达来进行模式匹配
• (2)features直接指定需要进行统计的基因名，优先于pattern

我们这里演示另外一个更底层的用法：

mt.genes=c('ENSSSCG00000018065','ENSSSCG00000018069',
     'ENSSSCG00000018075','ENSSSCG00000018078','ENSSSCG00000018080',
     'ENSSSCG00000018081','ENSSSCG00000018082','ENSSSCG00000018084',
     'ENSSSCG00000018086','ENSSSCG00000018087',
     'ENSSSCG00000018091','ENSSSCG00000018092','ENSSSCG00000018094')
head(rownames(sce.all))
kp=mt.genes %in% rownames(sce.all)
table(kp)
mt.genes=mt.genes[kp]
C<-GetAssayData(object = sce.all, slot = "counts")
percent.mito <- Matrix::colSums(C[mt.genes,])/Matrix::colSums(C)*100
sce.all <- AddMetaData(sce.all, percent.mito, col.name = "percent.mito")
sce.all[["percent.mito"]]
fivenum(sce.all@meta.data$percent.mito)

有意思的是，猪这个物种，它其实参考基因组的样式非常多，前面的ensembl数据库ID里面，是SSC的简称，有一些是CAF，这方面知识点应该是专门的做猪方面科研的小伙伴才能整理了。