植物的单细胞数据如何过滤线粒体基因

曾老师有一篇文章《猪的单细胞分析如何过滤线粒体基因》，其中介绍了猪的单细胞数据分析应该如何过滤线粒体基因，本期我们参考此文章来看看植物的单细胞数据分析如何过滤线粒体基因。

本期的数据来自《Single-Cell RNA Sequencing Resolves Molecular Relationships Among Individual Plant Cells》，据说是植物领域的第一篇单细胞转录组文章。

由于是植物的单细胞数据，不能像做人单细胞数据分析那样pattern = "^MT-"来去除线粒体的影响，起因是拟南芥的基因名没有特定的标记，因此我们得自己寻找基因列表。

下载线粒体基因列表

注释文件下载

首先我们进入EnsemblPlants，选中拟南芥的gff注释文件。

使用wget将注释文件下载到服务器上并解压。

wget http://ftp.ensemblgenomes.org/pub/plants/release-52/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.52.chromosome.Mt.gff3.gz
gunzip Arabidopsis_thaliana.TAIR10.52.chromosome.Mt.gff3.gz

将注释文件和单细胞数据基因名对应

> head(rownames(Seurat_object))
[1] "NAC001"    "ARV1"      "NGA3"      "DCL1"      "PPA1"      "AT1G01070"

通过查看我们发现Seurat对象中的基因名为symbol，因此在gff中我们也去找对应的ID类型。

在R中检测一下是否存在这个基因，果然是存在的。

> table(rownames(Seurat_object) %in% "CCB206")

FALSE  TRUE 
21964     1 

没有匹配上不代表这个注释文件有问题，一定要多试一些基因。

从注释文件中提取基因列表

接下来只需要提取我们需要的信息就好。

cat Arabidopsis_thaliana.TAIR10.52.chromosome.Mt.gff3 | grep -v "#" | awk -F "[\t=:;]" 'BEGIN{OFS="\t"}$3=="gene"{print $13}' > Ara_MTgenes.txt

• -F "[\t=:;]" 以这四个符号为分隔符
• 'BEGIN{OFS="\t"}13}' 只需要第三列为gene的行（忽略gene的大小写），并打印第十三列

打印多少列需要根据自己的注释文件进行修改，我需要的symbolID刚好在第13列。

wc查看共122个基因。

列表处理

我们在这里需要去除掉不包含在Seurat对象行名中的基因名。

MTgenes = unlist(read.csv(file = "Ara_MTgenes.txt",header = F))
table(MTgenes %in% rownames(Seurat_object))
MTgenes = MTgenes[MTgenes %in% rownames(Seurat_object)]

此处需要注意如果不使用unlist()函数，读取的不是基因列表，而是一个数据框，这会影响后续的%in%，而unlist()的作用就是将list数据变成字符串向量或者数字向量的形式。

过滤线粒体基因

添加线粒体信息

Seurat_object[["percent.mt"]] <- PercentageFeatureSet(
  Seurat_object, 
  features = MTgenes,
  )

可视化查看数据情况

如果有一些油滴里线粒体比例很高，而转录本很少，那可能是细胞已经破裂。

p1 <- FeatureScatter(Seurat_object, 
                     feature1 = "nCount_RNA", 
                     feature2 = "nFeature_RNA",
                     group.by = "sample")

p2 <- FeatureScatter(Seurat_object, 
                     feature1 = "nCount_RNA", 
                     feature2 = "percent.mt",
                     group.by = "sample")
p1 + p2 

过滤

# 过滤条件需要自行修改
Seurat_object <- subset(Seurat_object,
               subset = nFeature_RNA > 200 &
                 nFeature_RNA < 2500 &
                 percent.mt < 10)

报错

在写本期推文的时候并不是一番风顺，如果我没有进行列表处理去除掉不包含在Seurat对象行名中的基因名，而是直接运行添加线粒体信息，会报以下错误。

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'colSums': invalid character indexing

在github上看到一个issue讲的是这个，但不适用于本期内容。

以上。