Github带有全套代码分享的文献复现2025

生信技能树

发布于 2025-02-10 19:04:06

14300

代码可运行

运行总次数：0

代码可运行

我们生信技能树在前面给大家整理过大量的带有全套代码在 github 上分享的文献，还有一个专门的共享表格：自带图表复现代码及数据集的单细胞及多组学文章-V1 https://docs.qq.com/sheet/DT0FHR2NmbkVEVUJh。但是这样的资源是在太多了，这些优秀的代码可能被你收藏之后就再也未曾启封过。对大多数初学者来说，这些资源你可能想学习但心有力而不足，有些作者的代码实在是太复杂啦！我们秉承着一切分享学习的态度，即将开启2025年的新专辑《Github带有全套代码分享的文献复现2025》，我们以学习文献中的思路，代码技巧为主，带领大家攻破一篇篇文献，将其运行到自己的科研课题中。

下面第一篇学习的文章非常好，作者的数据，代码整理的非常详细见以前的分享：

代码分享|| 见过将代码整理成wiki资源的吗？这是篇极好的学习单细胞与scATAC-Seq组学以及联合分析的文献资源！

文献信息再次提供：

标题：A multi-omic single-cell landscape of human gynecologic malignancies
发表时间：2021年 12月发表在 Mol Cell （IF=14.5）杂志上
DOI：10.1016/j.molcel.2021.10.013
链接：https://linkinghub.elsevier.com/retrieve/pii/S1097-2765(21)00842-X
代码 github上也有：https://github.com/RegnerM2015/scENDO_scOVAR_2020/tree/v1.0.1
代码 wiki上的：https://github.com/RegnerM2015/scENDO_scOVAR_2020/wiki

这篇文章的分析也很有意思：

1、数据质控QC部分

作者是对每个样本单独进行了分析，之后对所有样本再进行合并。过滤的指标不是我们常见的过滤特定基因数或者count数的细胞，而是采用了一个MAD值，双细胞的判断则是用了两款软件 DoubletDecon 与 DoubletFinder，被同时鉴定为双细胞才是双细胞然后被过滤。

2、数据标准化部分

作者这里给出了不进行 UMI count 或者线粒体基因回归的原因，PC 个数的选择采用 JackStraw() 方法：

3、inferCNV 分析

这里作者对于正常参考 ref 的选择也不走寻常路：用了来自两个地方（TIMATE immune gene signature 和 PanglaoDB）的 signature 进行打分，在没有进行注释之前就找到免疫细胞打分高的cluster 当做正常对照，来进行拷贝数分析：

数据背景

11位未经治疗的患者接受了减瘤手术，肿瘤类型为子宫内膜或卵巢中发现的肿瘤。下表为每个病人的详细信息，每个病人同时做了单细胞转录组测序和单细胞 atac-seq，数据上传到了 GEO 中：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE173682：

下载地址：https://ftp.ncbi.nlm.nih.gov/geo/series/GSE173nnn/GSE173682/suppl/GSE173682_RAW.tar

scATAC-seq数据：

GSM5276944_3533EL_ATAC_fragments.tsv.gz 995.8 Mb
GSM5276945_3571DL_ATAC_fragments.tsv.gz 2.1 Gb
GSM5276946_36186L_ATAC_fragments.tsv.gz 1.6 Gb
GSM5276947_36639L_ATAC_fragments.tsv.gz 1.6 Gb
GSM5276948_366C5L_ATAC_fragments.tsv.gz 1.5 Gb
GSM5276949_37EACL_ATAC_fragments.tsv.gz 1.2 Gb
GSM5276950_38FE7L_ATAC_fragments.tsv.gz 535.7 Mb
GSM5276951_3BAE2L_ATAC_fragments.tsv.gz 1.5 Gb
GSM5276952_3E5CFL_ATAC_fragments.tsv.gz 1.4 Gb
GSM5276953_3CCF1L_ATAC_fragments.tsv.gz 1.3 Gb
GSM5276954_3E4D1L_ATAC_fragments.tsv.gz 1.4 Gb

scRNA-seq 数据：

GSM5276933_barcodes-3533EL.tsv.gz 29.2 Kb
GSM5276933_features-3533EL.tsv.gz 297.6 Kb
GSM5276933_matrix-3533EL.mtx.gz 44.3 Mb
GSM5276934_barcodes-3571DL.tsv.gz 40.1 Kb
GSM5276934_features-3571DL.tsv.gz 297.6 Kb
GSM5276934_matrix-3571DL.mtx.gz 60.1 Mb
GSM5276935_barcodes-36186L.tsv.gz 31.0 Kb
GSM5276935_features-36186L.tsv.gz 297.6 Kb
GSM5276935_matrix-36186L.mtx.gz 63.4 Mb
GSM5276936_barcodes-36639L.tsv.gz 40.8 Kb
GSM5276936_features-36639L.tsv.gz 297.6 Kb
GSM5276936_matrix-36639L.mtx.gz 60.1 Mb
GSM5276937_barcodes-366C5L.tsv.gz 35.3 Kb
GSM5276937_features-366C5L.tsv.gz 297.6 Kb
GSM5276937_matrix-366C5L.mtx.gz 55.1 Mb
GSM5276938_barcodes-37EACL.tsv.gz 40.3 Kb
GSM5276938_features-37EACL.tsv.gz 297.6 Kb
GSM5276938_matrix-37EACL.mtx.gz 67.3 Mb
GSM5276939_barcodes-38FE7L.tsv.gz 41.6 Kb
GSM5276939_features-38FE7L.tsv.gz 297.6 Kb
GSM5276939_matrix-38FE7L.mtx.gz 71.3 Mb
GSM5276940_barcodes-3BAE2L.tsv.gz 41.2 Kb
GSM5276940_features-3BAE2L.tsv.gz 297.6 Kb
GSM5276940_matrix-3BAE2L.mtx.gz 65.4 Mb
GSM5276941_barcodes-3CCF1L.tsv.gz 44.8 Kb
GSM5276941_features-3CCF1L.tsv.gz 297.6 Kb
GSM5276941_matrix-3CCF1L.mtx.gz 83.2 Mb
GSM5276942_barcodes-3E4D1L.tsv.gz 50.1 Kb
GSM5276942_features-3E4D1L.tsv.gz 297.6 Kb
GSM5276942_matrix-3E4D1L.mtx.gz 78.4 Mb
GSM5276943_barcodes-3E5CFL.tsv.gz 35.2 Kb
GSM5276943_features-3E5CFL.tsv.gz 297.6 Kb
GSM5276943_matrix-3E5CFL.mtx.gz 74.6 Mb

数据分析流程图

大致思路如下：

作者这里的策略是每个样本单独质控分析，然后再进行整合，可能主要的原因在于它使用了 Doublet 包进行双细胞去除！今天来学习一个样本的预处理！

单细胞部分：单个样本预处理

下面先来学习一个样本的处理看看，学会了一个样本的处理，其他的样本就都是类似了。这个样本预处理包括下面四个小步骤：

Preprocessing & QC
Feature selection, dimension reduction and clustering
inferCNV
SingleR reference-based annotation of cell types

0、需要提前安装的R包

当然，分析的时候你肯定需要做好环境配置的准备。如果下面这些包你安装有问题，然后自己又解决不了，那么后面的代码大概率也不适合你现在的R水平，可以提前放弃，或者R重新去练练，也可以考虑生信技能树的培训班，见文章末尾的友情宣传。

## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
options("repos"=c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"))
rm(list=ls())
library(scater)
library(dplyr)
library(Seurat)
library(patchwork)
library(SingleCellExperiment)
library(ComplexHeatmap)
library(ConsensusClusterPlus)
library(msigdbr)
library(fgsea)
library(tibble)
library(Signac)
library(ggplot2)
library(stringr)
library(SingleR)
library(devtools)
# DoubletDecon
# devtools::install_github('EDePasquale/DoubletDecon')
library(DoubletDecon)
library(DoubletFinder)
# BiocManager::install("EnsDb.Hsapiens.v86")
library(EnsDb.Hsapiens.v86)
library(xlsx)
library(GEOquery)
library(infercnv)

# rowr 这个包被cran删除了：https://cran.r-project.org/web/packages/rowr/index.html
# 可以下载以前的版本到本地安装
# wget https://cran.r-project.org/src/contrib/Archive/rowr/rowr_1.1.3.tar.gz
# bash 终端安装 -l 参数需要替换成自己的
# R CMD INSTALL -l /usr/local/software/miniconda3/envs/R4.4/lib/R/library rowr_1.1.3.tar.gz
library(rowr)

当然，还有一个比较重要的是目录的管理，你可能需要养成一个好习惯，整理好目录，我的目录如下，希望可以给你参考：

2021-GSE173682-单细胞多组学-卵巢和子宫内膜肿瘤
.
├── 00-文献资料
│   ├── 2021-GSE173682-卵巢和子宫内膜肿瘤-附件.pdf
│   ├── 2021-GSE173682-卵巢和子宫内膜肿瘤.pdf
│   ├── 2022-OMIX001928-同胞RNA+ATAC研究前列腺癌的谱系转变.pdf
│   ├── 2022-OV_And_Stemness-挖掘GSE173682.pdf
│   ├── CancerCell_同胞RNA+ATAC研究前列腺癌的谱系转变.pdf
│   └── GSE173682-GEO_Accession_viewer.pdf
├── 01-data        # 存放下载后的数据
├── 02-scRNA_Res   # 存放单细胞的分析结果
├── GSE173682      # 存放 GSE173682 的相关数据

到这里，你需要准备好环境，下载好数据和文献，代码，提前熟悉文献内容！

今天的学习目标：数据的QC部分

首先，准备好样本表型信息

# 样本表型信息
library(GEOquery)
gse <- "GSE173682"
gset <- getGEO(gse, destdir = './GSE173682/', getGPL = F)
a <- gset[[1]]
## 挑选一些感兴趣的临床表型。
pd <- pData(a)
colnames(pd)
# Patient_1_3533EL_RNA GSM5276933

这里我们先选择 Patient_1 即 GSM5276933 样本看看，作者的代码在：scENDO_scOVAR_2020-main\scRNA-seq Processing Scripts\Individual_Samples/Patient1_scRNA-seq.R中。

下载GEO数据

我这里在bash终端进行操作，下载到2021-GSE173682-单细胞多组学-卵巢和子宫内膜肿瘤/GSE173682

# 下载
wget -c https://ftp.ncbi.nlm.nih.gov/geo/series/GSE173nnn/GSE173682/suppl/GSE173682_RAW.tar
# 解压
tar -xvf GSE173682_RAW.tar

# 分别创建scATAC 和scRNA的目录
mkdir scATAC-Seq scRNA-Seq

# 将 scATAC文件放到 scATAC-Seq文件夹中
mv *fragments* scATAC-Seq/
ls

# 将 scRNA文件放到 scRNA-Seq文件夹中
mv *gz scRNA-Seq/
ls

现在开始整理 scRNA-Seq/ 目录中的文件为10X可读取的标准的三个文件：

cd scRNA-Seq/

## 创建每个样本文件夹
# 看一眼命令
ls *barcodes* |perl -ne 'chomp;/(.*)_(.*.gz)/;print"mkdir $1\n";'
# 运行
ls *barcodes* |perl -ne 'chomp;/(.*)_(.*.gz)/;print"mkdir $1\n";' |sh

# mv文件
ls *.gz |perl -ne 'chomp;/(.*)_(.*)-.*L(.*gz)/;print"mv $_ $1/$2$3\n";'
# 运行
ls *.gz |perl -ne 'chomp;/(.*)_(.*)-.*L(.*gz)/;print"mv $_ $1/$2$3\n";' |sh

整理成功后的目录如下：

你如果不是在服务器，那么上面的步骤你可以手动整理，或者见：人类妇科恶性肿瘤的多组学单细胞景观

采用MAD的细胞过滤

文章中的描述如下：每个样本单独进行质控，过滤细胞：

每个细胞中log后的UMI count数小于 2倍的 log后的UMI count数的MAD值，就会被过滤
每个细胞中表达的基因数：过滤方法同上
每个细胞中表达的线粒体基因比例：上面的小于改成大于

To enrich for high quality cells in each patient dataset, QC and doublet removal were performed for each patient dataset individually. First, outlier cells were defined in each of the following metrics: log(UMI counts) (> 2 MADs, low end), log(number of genes expressed) (> 2 MADs, low end) and log(percent mitochondrial read count +1) (> 2 MADs, high end) (McCarthyet al., 2017).

上面说的有点绕，代码运行来看看很快就会明白：

# Load the RNA dataset
# counts <- Read10X_h5(filename = "./filtered_feature_bc_matrix.h5")
counts <- Read10X("GSE173682/scRNA-Seq/GSM5276933/", gene.column = 2)

# Initialize the Seurat object with the raw (non-normalized data).
# genes not present in at least 3 cells are removed 
rna_raw <- CreateSeuratObject(counts = counts, min.cells = 3)
rna_raw

# 过滤前细胞数
PreQCNumCells <- length(colnames(rna_raw))
PreQCNumCells

# 计算线粒体基因表达比例，# 可能是13个线粒体基因
mito_genes <- grep("^MT-", rownames(rna_raw),ignore.case = T, value = T) 
print(mito_genes)
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
rna_raw[["percent.mt"]] <- PercentageFeatureSet(rna_raw, pattern = "^MT-")


# QC metrics: nCount_RNA, nFeature_RNA, and percent mitochrondial counts
# Outliers are >2 MADs
# Outliers函数：根据中位数绝对偏差（MAD）确定数值向量中的哪些值是异常值
# log(rna@meta.data$nCount_RNA)：对数转换用于稳定方差、减少偏态分布和突出显示数据中的小差异。
# lower" 表示函数将检测低于中位数的异常值（即下界异常值）
# nmads：异常值所需的最小MADs倍数
# nmads = 2：设置为2意味着如果一个值比中位数低2个MAD，它将被认为是一个异常值。
rna_raw@meta.data$nCount_RNA_outlier_2mad <- isOutlier(log(rna_raw@meta.data$nCount_RNA),log = F,type = "lower",nmads = 2)
table(rna_raw@meta.data$nCount_RNA_outlier_2mad)
rna_raw@meta.data$nFeature_RNA_outlier_2mad <- isOutlier(log(rna_raw@meta.data$nFeature_RNA),log = F,type = "lower",nmads = 2)
table(rna_raw@meta.data$nFeature_RNA_outlier_2mad)
rna_raw@meta.data$percent_mt_outlier_2mad <- isOutlier(log1p(rna_raw@meta.data$percent.mt),log = F,type = "higher",nmads = 2)
table(rna_raw@meta.data$percent_mt_outlier_2mad)

rna <- subset(rna_raw, subset = nCount_RNA_outlier_2mad == "FALSE" & 
                nFeature_RNA_outlier_2mad == 'FALSE' & 
                percent_mt_outlier_2mad == "FALSE")

# 过滤后的细胞数
PostQCNumCells <- length(colnames(rna))
PostQCNumCells
head(rna@meta.data)

过滤前 5697个细胞，过滤后5421个细胞，过滤掉了 276个细胞。

来看看小提琴分布：

p1 <- VlnPlot(rna_raw, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"), group.by = "orig.ident")
p2 <- VlnPlot(rna, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"), group.by = "orig.ident")
p1/p2

这种过滤方式实在是太温和了对吧！不要担心，作者后面还有双细胞过滤环节！