麻豆视频在播放,欧美图色日韩

文章標題

Single-cell RNA sequencing reveals distinct tumor microenvironmental patterns in lung adenocarcinoma

這篇文章通過單細胞RNA測序技術(shù)，對十例人類肺腺癌和十例正常組織進行了深入研究，發(fā)現(xiàn)肺癌細胞的轉(zhuǎn)錄組存在異質(zhì)性，反映了腫瘤的病理分級和基因突變途徑。同時，文章還發(fā)現(xiàn)存在兩種不同的微環(huán)境模式，即免疫激活的CP2E模式和惰性的N3MC模式。這兩種模式在預后方面具有不同的影響。此外，文章還發(fā)現(xiàn)微環(huán)境標記基因和簽名在整體腫瘤中具有預后價值。簡單來說，這篇文章通過單細胞RNA測序技術(shù)，提供了基于微環(huán)境的額外預后信息，并可能有助于預測治療反應和揭示未來治療方法的潛在目標細胞群。

雖然原文提供了代碼，但是作者之前用到的很多生信包和軟件已經(jīng)過時很久了，包括Seurat都已經(jīng)更新到V5了，所以今天就用自己的方法來復現(xiàn)一下文章中的fig.1中的圖。

Fig.1

一、準備安裝包和設(shè)置一些畫圖參數(shù)

library(Seurat)

library(dplyr)

library(reticulate)

library(sctransform)

library(cowplot)

library(ggplot2)

library(viridis)

library(tidyr)

library(magrittr)

library(reshape2)

library(readxl)

library(readr)

library(stringr)

# 后續(xù)繪制的 ggplot2 圖形都應用 cowplot 主題

theme_set(theme_cowplot())

#color scheme

use_colors <- c(Tumor ="brown2",# 腫瘤的顏色為棕色系

Normal ="deepskyblue2",# 正常的顏色為深天藍色系

G1 ="#46ACC8",# 細胞周期G1期的顏色為亮藍色

G2M ="#E58601",# 細胞周期G2/M期的顏色為橙色

S ="#B40F20",# 細胞周期S期的顏色為深紅色

Epithelial ="seagreen",# 上皮細胞的顏色為海洋綠

Immune ="darkgoldenrod2",# 免疫細胞的顏色為深金黃色

Stromal ="steelblue",# 間質(zhì)細胞的顏色為鋼藍色

p018 ="#E2D200",# p018標簽的顏色為檸檬黃

p019 ="#46ACC8",# p019標簽的顏色為亮藍色

p023 ="#E58601",# p023標簽的顏色為橙色

p024 ="#B40F20",# p024標簽的顏色為深紅色

p027 ="#0B775E",# p027標簽的顏色為深綠色

p028 ="#E1BD6D",# p028標簽的顏色為淺棕色

p029 ="#35274A",# p029標簽的顏色為深紫色

p030 ="#F2300F",# p030標簽的顏色為鮮紅色

p031 ="#7294D4",# p031標簽的顏色為淺藍色

p032 ="#5B1A18",# p032標簽的顏色為深棕色

p033 ="#9C964A",# p033標簽的顏色為淡黃色

p034 ="#FD6467"# p034標簽的顏色為粉紅色)

二、數(shù)據(jù)加載和質(zhì)控（到這之前基本都是用的原文代碼）

setwd("~/Data1/德國單細胞文章復現(xiàn)/01_初步處理/")

samples <- read_excel("../data/metadata/patients_metadata.xlsx", range = cell_cols("A:A")) %>% .$sample_id

for(iinseq_along(samples)){

# 讀取數(shù)據(jù)

assign(paste0("scs_data", i), Read10X(data.dir = paste0("../data/cellranger/", samples[i],"/filtered_feature_bc_matrix")))

}

### create seurat objects from cellranger files

# 使用for循環(huán)遍歷樣本列表中的每個樣本

for(iinseq_along(samples)) {

assign(paste0("seu_obj", i), CreateSeuratObject(

counts = eval(parse(text = paste0("scs_data", i))),# 將表達矩陣數(shù)據(jù)作為counts參數(shù)傳入，使用eval和parse函數(shù)動態(tài)獲取對應的變量

project = samples[i],# 將項目名稱作為project參數(shù)傳入

min.cells =3# 設(shè)置最小細胞數(shù)閾值為3，用于過濾低質(zhì)量細胞??

))

}

### merge data sets 合并多個Seurat對象為一個新的Seurat對象

seu_obj <- merge(seu_obj1,# 第一個原始Seurat對象

y = c(seu_obj2, seu_obj3, seu_obj4, seu_obj5, seu_obj6, seu_obj7, seu_obj8, seu_obj9, seu_obj10, seu_obj11, seu_obj12, seu_obj13, seu_obj14, seu_obj15, seu_obj16, seu_obj17, seu_obj18, seu_obj19, seu_obj20),add.cell.ids = samples,? project ="lung")

# add.cell.ids參數(shù)添加樣本標識，最后設(shè)置項目名稱為"lung"

seu_obj <- JoinLayers(seu_obj)

### calculate mitochondrial, hemoglobin and ribosomal gene counts

# 計算并添加MT基因組比例到Seurat對象中

seu_obj <- PercentageFeatureSet(seu_obj, pattern ="^MT-", col.name ="pMT")

# 計算并添加HBA/HBB基因組比例到Seurat對象中

seu_obj <- PercentageFeatureSet(seu_obj, pattern ="^HBA|^HBB", col.name ="pHB")

# 計算并添加RPS/RPL基因組比例到Seurat對象中,核糖體蛋白基因

seu_obj <- PercentageFeatureSet(seu_obj, pattern ="^RPS|^RPL", col.name ="pRP")

#### Data Filtering ####

## Before filtering

seu_obj_unfiltered <- seu_obj

qc_std_plot(seu_obj_unfiltered)

####After filtering####

seu_obj_unfiltered <- readRDS("seurat_objects/all_unfiltered.RDS")

#過濾

seu_obj <- subset(seu_obj_unfiltered, subset = nFeature_RNA > nFeature_lower & nFeature_RNA < nFeature_upper & nCount_RNA > nCount_lower & nCount_RNA < nCount_upper & pMT < pMT_upper & pHB < pHB_upper)

SuppFig1A,Before filtering

SuppFig1B

三、 Data normalization/harmony去批次 /降維聚類

seu_obj<- NormalizeData(seu_obj, normalization.method="LogNormalize",scale.factor=1e4)??

seu_obj<- FindVariableFeatures(seu_obj)?

#?FindVariableFeatures

seu_obj<- ScaleData(seu_obj)

#?ScaleData

seu_obj<- RunPCA(seu_obj)

#?RunPCA（挑選顯著PC，最好run一下，避免引入噪音）

ElbowPlot(seu_obj,ndims=50) ####?

#Seurat 5版本有斷層數(shù)據(jù)需要連接下

# seu_obj <- JoinLayers(seu_obj)

library(harmony)??

###Harmony

options(repr.plot.height=2.5, repr.plot.width = 6)

seu_obj_harmony<- seu_obj %>%

RunHarmony(group.by.vars="orig.ident", plot_convergence = TRUE)??

#### Dimensionality reduction降維聚類####

seu_obj_harmony<- FindNeighbors(seu_obj_harmony, reduction = "harmony",dims=1:15)?

seu_obj_harmony<- FindClusters(seu_obj_harmony, resolution = 0.2)

seu_obj_harmony<- RunUMAP(seu_obj_harmony,? dims = 1:15, reduction="harmony")

print(DimPlot(seu_obj_harmony, reduction = "umap",label = T) + labs(title = paste0("resolution: ", 0.2)))

harmony去批次后的降維聚類圖是不是和文獻中的UMAP圖形狀大致吻合

四、Main cell type annotation/大類細胞注釋

mainmarkers<c("PECAM1","VWF","ACTA2","JCHAIN","MS4A1","PTPRC","CD68","KIT","EPCAM","CDH1","KRT7","KRT19")

DotPlot(seu_obj_harmony, features = mainmarkers, group.by ="RNA_snn_res.0.2") +??

coord_flip() +??

?scale_color_viridis()

ggsave2("DotPlot_mainmarkers.png", path ="output2/annotation", width =30, height =8, units ="cm")

DimPlot(seu_obj_harmony, group.by ="RNA_snn_res.0.2", label = T, label.size =5)

ggsave2("DimPlot_all_clusters.png", path ="output/annotation", width =20, height =20, units ="cm")

Idents(seu_obj_harmony) <- seu_obj_harmony$RNA_snn_res.0.2

#這里三個分成大類免疫、上皮、基質(zhì)

# saveRDS(seu_obj, file = "seurat_objects/SCT_snn_res.0.085.RDS")

# seu_obj<- readRDS("seurat_objects/SCT_snn_res.0.085.RDS")

annotation_curated_main <-read_excel("../data/curated_annotation/curated_annotation_main改_harmony.xlsx")

new_ids_main <- annotation_curated_main$main_cell_type

names(new_ids_main) <- levels(seu_obj_harmony)

seu_obj_harmony <- RenameIdents(seu_obj_harmony, new_ids_main)

seu_obj_harmony@meta.data$main_cell_type<- Idents(seu_obj_harmony)

DimPlot(seu_obj_harmony, group.by ="main_cell_type")

DimPlot(seu_obj_harmony, group.by ="RNA_snn_res.0.2",label = T)

# Add metadata

metatable <- read_excel("../data/metadata/patients_metadata.xlsx")

metadata <- FetchData(seu_obj_harmony,"orig.ident")

metadata$cell_id<- rownames(metadata)

metadata$sample_id<- metadata$orig.ident

metadata <- left_join(x = metadata, y = metatable, by ="sample_id")

rownames(metadata) <- metadata$cell_id

seu_obj_harmony <- AddMetaData(seu_obj_harmony, metadata = metadata)

我用harmony去批次整合的圖

文獻fig.1D的原圖（肉眼可以看到細胞形狀和細胞數(shù)量大致相同）

五、Cell cycle scoring/細胞周期評分

### add cell cycle, cc.genes loaded with Seurat

s.genes <- cc.genes$s.genes

g2m.genes <- cc.genes$g2m.genes

score_cc <-function(seu_obj_harmony){??

seu_obj_harmony <- CellCycleScoring(seu_obj_harmony, s.genes, g2m.genes)??

seu_obj_harmony@meta.data$CC.Diff <- seu_obj_harmony@meta.data$S.Score - seu_obj_harmony@meta.data$G2M.Score

return(seu_obj_harmony)

}

seu_obj_harmony <- score_cc(seu_obj_harmony)

FeatureScatter(seu_obj_harmony,"G2M.Score","S.Score", group.by ="Phase") +

coord_fixed(ratio =1)

可以看到大部分細胞周期評分都很低可以不用去除Cell cycle相關(guān)基因

六、 plots for figure 1

DimPlot(seu_obj_harmony,group.by="tissue_type", cols = use_colors)

ggsave2("Fig1B.png",?path?="../results",?width?=15,?height?=15,?units?="cm")

DimPlot(seu_obj_harmony,group.by="patient_id", cols = use_colors, pt.size =1)

ggsave2("Fig1C.png",?path?="../results",?width?=15,?height?=15,?units?="cm")

DimPlot(seu_obj_harmony,group.by="main_cell_type", cols = use_colors, pt.size =1)

ggsave2("Fig1D_umap.png",?path?="../results",?width?=15,?height?=15,?units?="cm")

cell_types <- FetchData(seu_obj_harmony, vars = c("sample_id","main_cell_type","tissue_type")) %>%

mutate(main_cell_type = factor(main_cell_type, levels = c("Stromal","Immune","Epithelial"))) %>%

mutate(sample_id = factor(sample_id, levels = rev(c("p018t","p019t","p023t","p024t","p027t","p028t","p030t","p031t","p032t","p033t","p034t","p018n","p019n","p027n","p028n","p029n","p030n","p031n","p032n","p033n","p034n"))))

ggplot(data = cell_types) +?

geom_bar(mapping = aes(x = sample_id, fill = main_cell_type, ), position ="fill", width =0.75) +?

scale_fill_manual(values = use_colors) +??

coord_flip()

ggsave2("Fig1D_barplot.pdf", path ="../results", width =15, height =30, units ="cm")

remove(seu_obj_harmony)