复现文章：Analyses of metastasis-associated genes in IDH wild-type glioma
数据来源：GSE84465
图1
image.png
数据筛选条件
image.png
image.png
除了min.cells min.features percent.mt，还有Pearson相关系数的指标（Fig S1）
根据Pearson相关系数大于0.4这个指标，去除了Tumor_BT_S1和Tumor_BT_S4两个样本，只针对剩余的6个样本分析
Fig S1
image.png
第一步：下载GSE84465_GBM_All_data.csv.gz和Metadata的SraRunTable.txt文件
第二步代码
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
a<-read.table("GSE84465_GBM_All_data.csv.gz")
head(rownames(a))
tail(rownames(a),10)
a=a[1:(nrow(a)-5),] #最后5行行名异常，剔除
b<-read.table("SraRunTable.txt",sep = ",",header = T) #样本信息
table(b$PATIENT_ID)
table(b$tissue)
table(b$tissue, b$PATIENT_ID)
new.b <- b[,c("Plate.ID","well","tissue","PATIENT_ID")]
new.b$sample <- paste0("X",b$Plate.ID,".",b$well)
head(new.b)
identical(colnames(a),new.b$sample)
pbmc <- CreateSeuratObject(counts = a)
head(pbmc@meta.data)
patient<-new.b$PATIENT_ID
pbmc <- AddMetaData(object = pbmc, metadata = patient, col.name = 'PATIENT_ID')
tissue<-new.b$tissue
pbmc <- AddMetaData(object = pbmc, metadata = tissue, col.name = 'Tissue')
sample<-paste0(tissue,"_",patient)
pbmc <- AddMetaData(object = pbmc, metadata = sample, col.name = 'Sample')
head(pbmc@meta.data)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
table(pbmc$PATIENT_ID,pbmc$Tissue)
###########################
#计算Pearson相关系数
pbmc.filt <- subset(pbmc, subset = PATIENT_ID=="BT_S6"&Tissue=="Periphery") #BT_S1 2 4 6 Tumor Periphery
dim(pbmc.filt)
plot1 <- FeatureScatter(pbmc.filt, feature1 = "nCount_RNA", feature2 = "percent.mt",group.by = "PATIENT_ID",pt.size=1.5)
plot2 <- FeatureScatter(pbmc.filt, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",group.by = "PATIENT_ID",pt.size=1.5)
plot1 + plot2
############################
selected_f <- rownames(pbmc)[Matrix::rowSums(pbmc@assays$RNA@counts > 0 ) > 3]
pbmc.filt <- subset(pbmc, subset = nFeature_RNA > 50 & nFeature_RNA < 6000 & percent.mt < 0.05 &
Sample%in%c("Periphery_BT_S1","Periphery_BT_S2","Periphery_BT_S4","Periphery_BT_S6","Tumor_BT_S2","Tumor_BT_S6"),
features = selected_f)
dim(pbmc.filt)
VlnPlot(object = pbmc.filt, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),group.by = "Sample")
image.png
image.png
#标准化
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
### 鉴定高变异基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 1500) #变异最大的1500个基因
#输出特征方差图
top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10,repel = TRUE)
plot1+plot2
image.png
图2
image.png
#标准化
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
### 鉴定高变异基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 500) #变异最大的1500个基因
#输出特征方差图
top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10,repel = TRUE)
plot1+plot2
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(object = pbmc, dims = 1:2, reduction = "pca",nfeatures = 20) #4个PC 20个基因
DimPlot(pbmc, reduction = "pca",group.by="Sample")
DimHeatmap(pbmc, dims = 1:2, cells = 500, balanced = TRUE,nfeatures = 30,ncol=2)
pbmc <- JackStraw(object = pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)
JackStrawPlot(object = pbmc, dims = 1:20,reduction = "pca")
ElbowPlot(pbmc,reduction="pca") #根据ElbowPlot确定PC数量
##TSNE聚类分析
pbmc <- FindNeighbors(pbmc, dims = 1:20)
pbmc <- FindClusters(pbmc, resolution = 0.5)
#tSNE降维
pbmc <- RunTSNE(object = pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "tsne",label = TRUE,pt.size = 1)
image.png
image.png
image.png
image.png
image.png
table(pbmc@meta.data$Tissue,pbmc@meta.data$seurat_clusters)
image.png
### 差异表达分析
# 找到每一个cluster当中的marker，并且只展示阳性的marker。
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
head(pbmc.markers, n = 10)
write.table(pbmc.markers,file = 'pbmc.markers.txt',sep = '\t',row.names=F,quote=F)
library("tidyverse")
sig.markers<- pbmc.markers %>% select(gene,everything()) %>%
subset(p_val_adj<0.05 & abs(pbmc.markers$avg_log2FC)>1)
dim(sig.markers)
write.table(sig.markers,file="sigmarkers.xls",sep="\t",row.names=F,quote=F)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
#差异基因可视化
VlnPlot(object = pbmc, features = c("EGFR", "CCL3","AGXT2L1","DCN","GPR17","MOG","SYT1"),group.by = "seurat_clusters",pt.size = 1)
image.png
Cluster 2 3 6被定义为metastatic tumor cell，比较metastatic和non-mentastatic之间的基因表达差异
metastatic.markers <- FindMarkers(pbmc, ident.1 = c(2,3,6), ident.2 = c(0,1,4,5,7,8,9,10,11,12), min.pct = 0.25)
head(metastatic.markers,10)
sig.metastatic.markers<- metastatic.markers %>% select(colnames(metastatic.markers),everything()) %>%
subset(p_val_adj<0.05 & abs(metastatic.markers$avg_log2FC)>1)
dim(sig.metastatic.markers)
write.table(sig.metastatic.markers,file="sigmetastaticmarkers.xls",sep="\t",row.names=F,quote=F)
网友评论