search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
124 stars 30 forks source link

scanpy和Seurat单细胞分析对比 #4111

Closed ixxmu closed 1 year ago

ixxmu commented 1 year ago

https://mp.weixin.qq.com/s/JtBn0iEW2N9GBK5Oj0UJqQ

ixxmu commented 1 year ago

scanpy和Seurat单细胞分析对比 by 生信技能树

曾老师已经是V4的版本的Seurat流程做出来了大量的公共数据集的单细胞转录组降维聚类分群流程，100多个公共单细胞数据集全部的处理，链接：https://pan.baidu.com/s/1MzfqW07P9ZqEA_URQ6rLbA?pwd=3heo

但是最近Seurat版本问题，导致好多人使用我们的标准代码失败了，首先是可以参考我们前面的降级方式：假如你不喜欢最新版的Seurat包的单细胞理念，其次也可以干脆就转为Python的单细胞scanpy流程

scanpy和seurat分析代码比较

尝试把曾老师的单细胞seurat分析的代码转换成scanpy版本的，包括样品读取，质控，harmony去除批次效应，降维聚类，marker鉴定。

seurat版本的R代码

先看看seurat版本的R代码，里面调用了另外几个文件，lib.R，qc.R，harmony.R，check-all-markers.R，这几个文件以前应该都有分享过.scRNA_scripts文件夹的r代码都在上面的百度云网盘链接（https://pan.baidu.com/s/1MzfqW07P9ZqEA_URQ6rLbA?pwd=3heo）里面哦：

rm(list=ls())
options(stringsAsFactors = F) 
source('scRNA_scripts/lib.R')

###### step1:导入数据 ######   
# 付费环节 800 元人民币

dir='GSE189357_RAW/output/' 
samples=list.files( dir  )
samples 

library(data.table)
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro) 
  sce=CreateSeuratObject(counts =  Read10X(file.path(dir,pro) ) ,
                         project =  pro ,
                         min.cells = 5,
                         min.features = 300,)
  
  return(sce)
})

names(sceList)  

library(stringr)
# samples=gsub('.txt.gz','',str_split(samples,'_',simplify = T)[,2])
samples

names(sceList) =  samples

sce.all <- merge(sceList[[1]], y= sceList[ -1 ] ,
                 add.cell.ids =  samples) 

as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all@meta.data$orig.ident) 
 

###### step2:QC质控 ######
dir.create("./1-QC")
setwd("./1-QC")
# 如果过滤的太狠，就需要去修改这个过滤代码
source('../scRNA_scripts/qc.R')
sce.all.filt = basic_qc(sce.all)
print(dim(sce.all))
print(dim(sce.all.filt))
setwd('../')


sp='human'

###### step3: harmony整合多个单细胞样品 ######
dir.create("2-harmony")
getwd()
setwd("2-harmony")
source('../scRNA_scripts/harmony.R')
# 默认 ScaleData 没有添加"nCount_RNA", "nFeature_RNA"
# 默认的
sce.all.int = run_harmony(sce.all.filt)
setwd('../')

###### step4: 降维聚类分群和看标记基因库 ######
# 原则上分辨率是需要自己肉眼判断，取决于个人经验
# 为了省力，我们直接看  0.1和0.8即可
table(Idents(sce.all.int))
table(sce.all.int$seurat_clusters)
table(sce.all.int$RNA_snn_res.0.1) 
table(sce.all.int$RNA_snn_res.0.8) 

getwd()
dir.create('check-by-0.1')
setwd('check-by-0.1')
sel.clust = "RNA_snn_res.0.1"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 

source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

dir.create('check-by-0.8')
setwd('check-by-0.8')
sel.clust = "RNA_snn_res.0.8"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

dir.create('check-by-0.5')
setwd('check-by-0.5')
sel.clust = "RNA_snn_res.0.5"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

last_markers_to_check

scanpy版本的python代码

安装各种库的时候，直接pip install scanpy即可

如果是在jupyter notebook或者jupyter lab里，需要在前面加上！，命令行里就不用

!pip install scanpy

#安装导入需要的模块
import scanpy as sc
import numpy as np 
import pandas as pd 
from glob import glob
import re
import seaborn as sns
import matplotlib.pyplot as plt 
from tqdm import tqdm
import os
from collections import Counter
import time
import cosg as cosg
%matplotlib inline
# 开始计时
a = time.perf_counter()

scanpy下的read_10x_mtx函数读取数据，传入文件夹路径，保证文件夹下还是这三个文件即可

循环读取9个单细胞数据，用字典进行存储，key为样本名，value为scanpy读取后的对象

用scanpy下的concat函数进行多个样本的合并

scanpy对象结构参考下图

adata.var是行为gene，列为统计指标的矩阵，进行完各种分析之后，大概如下图，刚读完数据的时候，只有行名，没有列

adata.obs是行为细胞，列为统计指标的矩阵adata.uns里面存放的是非结构化的数据，都是字典的形式

比如sample color，存放的是每个样品对应的颜色，pca里面存放有降维后的主成分的信息，louvain 0.01_colors里面存放的就是0.01的分辨率下每个cluster对应的颜色

adata.X里面存放的是表达矩阵，行为cell，列为gene

### Step1 批量读取单细胞数据
# 字符串前面加上r之后，路径里的/和\就不用特意修改了
path = r'F:\新建文件夹\2022-GSE189357-LUAD-单细胞-疾病进展\GSE189357_RAW\output'
files = glob(path + '\*')

# tqdm用来显示进度条，不加也可以
data_dic = {}
for i in tqdm(files):
    sample_name = i.split('\\')[-1]
    # cache=True 会存储读取文件时候的缓存，下次再读取就很快了
    adata = sc.read_10x_mtx(i,var_names='gene_symbols',cache=True)
    adata.var_names_make_unique()
    adata.obs_names_make_unique()
    data_dic[sample_name] = adata
    
adata = sc.concat(data_dic,label='sample',axis=0)
adata.obs_names_make_unique()

os.getcwd()

work_dir = r'F:新建文件夹/'
os.chdir(work_dir)

qc部分，参考曾老师seurat里qc.R脚本写的

这里只考虑了human gene，都是大写的

#计算比例和标准质控  
def basic_qc(adata):
    # 计算线粒体基因比例
    adata.var['mt'] = adata.var_names.str.startswith('MT-')
    sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
    mt_gene = adata.var.index[adata.var_names.str.startswith('MT-')]
    print('线粒体gene：',mt_gene)
    # 计算核糖体基因比例
    adata.var['rp'] = [True if i.startswith('RPL') or i.startswith('RPS') else False for i in adata.var_names]
    sc.pp.calculate_qc_metrics(adata, qc_vars=['rp'], percent_top=None, log1p=False, inplace=True)
    rp_gene = [i for i in adata.var_names if i.startswith('RPL') or i.startswith('RPS')]
    print('核糖体gene：',rp_gene)
    # 计算红血细胞基因比例
    adata.var['hb'] = [True if i.startswith('HB') and not i.startswith('HBP') else False for i in adata.var_names]
    sc.pp.calculate_qc_metrics(adata, qc_vars=['hb'], percent_top=None, log1p=False, inplace=True)
    hb_gene = [i for i in adata.var_names if i.startswith('HB') and not i.startswith('HBP')]
    print('红血细胞gene：',hb_gene)
    
    # 可视化细胞各gene比例
    if not os.path.isdir('qc'):
        os.mkdir('qc')
    
    sc.pl.violin(adata,keys=['n_genes_by_counts','total_counts'],groupby='sample',show = False)
    plt.savefig('qc/vlnplot1.pdf')
    sc.pl.violin(adata,keys=['pct_counts_mt','pct_counts_rp','pct_counts_hb'],groupby='sample',show = False)
    plt.savefig('qc/vlnplot2.pdf')
    sc.pl.scatter(adata,x='total_counts', y='n_genes_by_counts',color='sample',show = False)
    plt.savefig('qc/scatterplot.pdf')
    sc.pl.highest_expr_genes(adata,n_top=20,show = False)
    plt.savefig('qc/TOP20_most_expressed_gene.pdf')
    
    # 过滤
    print('过滤前：',adata.shape)
    adata = adata[adata.obs.n_genes_by_counts < 6000, :]
    adata = adata[adata.obs.pct_counts_mt < 20, :]
    adata = adata[adata.obs.pct_counts_rp > 3, :]
    adata = adata[adata.obs.pct_counts_hb < 1, :]
    print('过滤后：',adata.shape)
    
    # 可视化过滤后的结果
    sc.pl.violin(adata,keys=['n_genes_by_counts','total_counts'],groupby='sample',show = False)
    plt.savefig('qc/vlnplot1_filtered.pdf')
    sc.pl.violin(adata,keys=['pct_counts_mt','pct_counts_rp','pct_counts_hb'],groupby='sample',show = False)
    plt.savefig('qc/vlnplot2_filtered.pdf')
    return adata
#harmony整合

def run_harmony(adata):
    # normalize
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    adata.raw = adata
    # 找高变gene
    sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
    # 只保留高变gene
    adata = adata[:, adata.var.highly_variable]
    # 根据线粒体gene进行regress
    sc.pp.regress_out(adata, 'pct_counts_mt')
    # scale
    sc.pp.scale(adata)
    # pca
    sc.tl.pca(adata, svd_solver='arpack')
    sc.external.pp.harmony_integrate(adata,key='sample')
    
    sc.pp.neighbors(adata, n_neighbors=10, n_pcs=20,use_rep='X_pca_harmony')
    sc.tl.umap(adata)
    
    # 设置不同分辨率
    for i in [0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1]:
        sc.tl.louvain(adata,resolution=i,key_added='louvain {}'.format(i))
    if not os.path.isdir('harmony'):
     os.mkdir('harmony')
    sc.pl.umap(adata,color=['louvain 0.01','louvain 0.1','louvain 0.2'],show = False)
    plt.savefig('harmony/resolution_low.pdf')
    sc.pl.umap(adata,color=['louvain 0.5','louvain 0.8','louvain 1'],show = False)
    plt.savefig('harmony/resolution_high.pdf')
    return adata
def run_harmony(adata):
    # normalize
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    adata.raw = adata
    # 找高变gene
    sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
    # 只保留高变gene
    adata = adata[:, adata.var.highly_variable]
    # 根据线粒体gene进行regress
    sc.pp.regress_out(adata, 'pct_counts_mt')
    # scale
    sc.pp.scale(adata)
    # pca
    sc.tl.pca(adata, svd_solver='arpack')
    sc.external.pp.harmony_integrate(adata,key='sample')
    
    sc.pp.neighbors(adata, n_neighbors=10, n_pcs=20,use_rep='X_pca_harmony')
    sc.tl.umap(adata)
    
    # 设置不同分辨率
    for i in [0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1]:
        sc.tl.louvain(adata,resolution=i,key_added='louvain {}'.format(i))
    if not os.path.isdir('harmony'):
     os.mkdir('harmony')
    sc.pl.umap(adata,color=['louvain 0.01','louvain 0.1','louvain 0.2'],show = False)
    plt.savefig('harmony/resolution_low.pdf')
    sc.pl.umap(adata,color=['louvain 0.5','louvain 0.8','louvain 1'],show = False)
    plt.savefig('harmony/resolution_high.pdf')
    return adata

结果展示

adata = basic_qc(adata)
adata = run_harmony(adata)
check_markers(0.1)
check_markers(0.3)
check_markers(0.8)

运行完耗时1021s

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：