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

当我们进行非模式物种的单细胞分析时，难免会遇到各种问题，其中，基因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([e-mail]$percent.mito)

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

文末友情推荐

做教学我们是认真的，如果你对我们的马拉松授课（直播一个月互动教学）有疑问，可以看完我们从2000多个提问互动交流里面精选的200个问答！ 2021第二期_生信入门班_微信群答疑整理，以及 2021第二期_数据挖掘班_微信群答疑笔记

与十万人一起学生信，你值得拥有下面的学习班：