新TCGA+文献复现里的几种算法

原创

用户10758803

发布于 2024-03-10 01:53:58

3000

差异分析的起点：count矩阵，只能用count数据做差异分析

代码和图片均来自生信技能树小洁老师

reads计数数据（测序的短片段），会匹配到基因。若匹配到，则匹配到的基因会count+1。（一个基因对应4个read，即count为4）

Gtex：正常样本的组织?

TCGA 正常组织样本少，可以与Gtex联合。在UCSC xena网站上

（点击UCSC Toil RNAseq Recompute Compendium）

差异分析

输入数据：exp、分组向量、proj（做文件前缀、项目名称)

如何检查自己的代码有无搞反？

x = rownames(DEG1)[1]
dat = log2(cpm(exp)+1)
boxplot(dat[x,]~Group)
DEG1[1,]

做差异分析时用count，差异分析可视化时用cpm、tpm

构建模型时，样本数据多多益善（>100个）

基因过滤后就开始用LogCPM或LogTPM数据

生存信息，去xena里找，xena里也有临床信息，但临床信息已经用TCGAlink下载过。以病人iid列连接在一起

表达矩阵与临床信息需要匹配，否则没办法把一个基因当作一个临床因素去处理

KM曲线

可以直观展示生存率和死亡率,有p值，展示组间生存率变化的比较

log_rank_test

log_rank_test：批量展示一群基因的p值，没有图，只有计算结果。结果为一组有名字的向量。

用于批量计算p值

library(survival)
library(survminer)
log_rank_p <- apply(exprSet , 1 , geneKM)
  #genekm是老师打包好的函数，计算出基因的p值
  log_rank_p=sort(log_rank_p)
  save(log_rank_p,file = logrankfile)

  tinyarray::surv_KM()#老师写的包里的函数

批量单因素cox

直接展示一个表格，有p值、HR（重要统计指标。大于1则表示危险因素）

不用于估计生存率，用于因素分析，找到某一个因素对结局事件发生的贡献度

只有他离散数据和连续数据都可以接受。

library(survival)
library(survminer)
cox_results <-apply(exprSet , 1 , genecox)#genecox也是老师写的函数
  cox_results=as.data.frame(t(cox_results))#转置并转化成数据框
  save(cox_results,file = coxfile)
  tinyarray::surv_cox()#老师写的包里的函数

做模型：挑出关键基因，得到公式，给每一个病人计算风险分数

lasson回归、COX多因素分析、随机森林、支持向量随机

缩小基因的数量，得到公式，得到风险分数

最要学会的数据整理方法：TCGA_2里的 5.sur.dat

GBM_ER里 IDH突变和MTMG启动子甲基化、1p19q共缺失的数据整理

GSEA和GSVA

https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/mxqi4z#5f199b96

怎么做？

共同的数据准备

1.数据准备：exp（一列一个样本，一行一个基因），Group（样本的分组），limma差异分析结果（exp差异分析的结果，只要log_FC即可）

2.数据包misgdbr（这里其实是构建一个文库，然后再把差异基因在这里面取子集？

KEGG_df = msigdbr(species = "Homo sapiens",category = "C2",subcategory = "CP:KEGG") %>% 
  dplyr::select(gs_exact_source,gene_symbol)
head(KEGG_df)
## # A tibble: 6 x 2
##   gs_exact_source gene_symbol
##   <chr>           <chr>      
## 1 hsa02010        ABCA1      
## 2 hsa02010        ABCA10     
## 3 hsa02010        ABCA12     
## 4 hsa02010        ABCA13     
## 5 hsa02010        ABCA2      
## 6 hsa02010        ABCA3

GO_df = msigdbr(species = "Homo sapiens",category = "C5") %>% 
  dplyr::select(gene_symbol,gs_exact_source,gs_subcat)
dim(GO_df)
## # A tibble: 6 x 2
##   gs_exact_source gene_symbol
##   <chr>           <chr>      
## 1 GO:0004645      GDPGP1     
## 2 GO:0004645      MTAP       
## 3 GO:0004645      PYGB       
## 4 GO:0004645      PYGL       
## 5 GO:0004645      PYGM       
## 6 GO:0004645      TYMP

GSEA(基因集富集分析)

GSEA需要数据：将logFC从小到大排序，组成一个向量，并且给每一个向量命上对应的基因名。构成的向量名叫ge，还有GO_df

ge = deg$logFC
names(ge) = rownames(deg)
ge = sort(ge,decreasing = T)
head(ge)
##     OLR1   COL1A1    OLFM4     H3C8   CRISP3  ZFP36L2 
## 3.835314 3.790745 3.427338 3.180043 2.956932 2.909744
em <- GSEA(ge, TERM2GENE = GO_df)
#画个图来看看
gseaplot2(em, geneSetID = 1, title = em$Description[1])

从下往上看：

Ranked list metric：柱状图叠加，所有基因的叠加，展示所有基因的logFC变化。横坐标代表的时logFC

颜色图：按照logFC从大到小颜色分配，负值代表蓝色，正的代表红色

条形码：一条竖线展示一个基因，和下面的logFC对应。空白的地方代表这个位置的基因没有属于这条通路的

running enchscore展示的曲线：每发现一条基因在这个通路上，runnscore会增加一些？曲线向下，下调基因富集？

GSVA（基因集差异分析）

评估不同代谢通路在不同样本里

其实就像将原来exp的行名转化为通路的名字。这里需要将exp的行名的基因富集到具体的通路上，再将通路换到exp行名上，生成一个新的表达矩阵kegg_ES/go_ES同理,接来来就是做limma差异分析，接下来就可以用热图、火山图等可视化

kegg_list = split(KEGG_df$gene_symbol,KEGG_df$gs_exact_source)#将每个通路具体的基因找出来
lapply(kegg_list[1:3], head)
## $hsa00010
## [1] "ACSS1" "ACSS2" "ADH1A" "ADH1B" "ADH1C" "ADH4" 
## 
## $hsa00020
## [1] "ACLY" "ACO1" "ACO2" "CS"   "DLAT" "DLD" 
## 
## $hsa00030
## [1] "ALDOA" "ALDOB" "ALDOC" "DERA"  "FBP1"  "FBP2"
KEGG_ES <- gsva(expr=exp, 
               gset.idx.list=kegg_list, 
               parallel.sz=32) #自己电脑parallel.sz写5就好，线程数
KEGG_ES[1:4,1:4]
##          GSM1366348 GSM1366349 GSM1366350  GSM1366351
## hsa00010 0.09619189 -0.2030580 -0.2575823  0.13014453
## hsa00020 0.41117392 -0.4726020 -0.4338775 -0.26959554
## hsa00030 0.03504000 -0.4179045 -0.4280300 -0.09500726
## hsa00040 0.16960156 -0.1391068 -0.2706679 -0.02557599
design = model.matrix(~Group)
fit = lmFit(KEGG_ES, design)
fit = eBayes(fit)
DEG = topTable(fit, coef = 2, number = Inf)
head(DEG)

随机森林

必须要有测试集和训练集

TCGA联合分析:根据某个基因是否突变分组比较某基因的表达量

TCGA多组学数据，桥架为样本ID，不管他是测了突变信息还是测了RNAseq。

根据病人ID就可以把同一个病人的不同数据联合在一起。

表达矩阵里有的样本在突变数据maf里不一定有，要把没有突变的病人去掉

得到的小提琴图表示：VHL基因的突变是否影响PBRM1的表达

如果影响则小提琴组会有较大的差别

分组：基因表达量的高级？是否甲基化？是否突变？

任意基因的相关性

可以将分组（正常样本和肿瘤样本）与基因的相关性联系

几种算法（免疫丰度、免疫细胞亚型

1.免疫亚型鉴定和可视化

https://www.yuque.com/xiaojiewanglezenmofenshen/bsgk2d/ab5gvv#6d0c41bc)

结合GBM_ER的文档3里的代码immune里来看

准备好exp？（取过log和未取过log的数据都可以，但是TPM，FPKM，cpm不行）

再用函数callEnsemble

res0 <- callEnsemble(X =exp, geneids = 'symbol')
res0

输出数据，这里的bestcall即是样本的算出的最佳的免疫亚型，后面的1-6应该是每种免疫亚型的得分?最高的得分即使bestcall？

##    SampleIDs BestCall            1            2            3            4
## 1        XY1        4 1.121173e-07 1.873900e-06 6.311234e-02 0.9027497470
## 2        XY2        3 4.243225e-02 2.309184e-06 5.495564e-01 0.0084770960
## 3        XY3        4 3.095071e-04 1.029907e-06 1.635270e-01 0.9063920975
## 4        XY4        4 1.154656e-04 3.823049e-07 2.888787e-02 0.8831390142


##               5            6
## 1  5.132385e-02 5.121503e-05
## 2  1.000665e-04 2.261875e-04
## 3  5.731729e-03 1.741630e-04
## 4  1.236814e-03 7.847964e-05

2.我们只需要1，2列。

3.（小洁老师的链接里：可以检验一下结果是否可靠，即geneMatchErrorReport告诉你有多少基因是在输入数据中不存在的，如果比例过高，那最终计算出来的亚型就不那么可靠）

4.作图了

library(dplyr)
dat2 <- dat %>%
 group_by(group, BestCall) %>% 
 summarise(count = n())

ggplot() + geom_bar(data =dat2, aes(x = group, y = count, fill = BestCall),
 stat = "identity",
 position = "fill")+
  theme_classic()+
  scale_fill_paletteer_d("RColorBrewer::Set2")

##2. 免疫细胞丰度

箱线图：文献代码：https://www.yuque.com/xiaojiewanglezenmofenshen/bsgk2d/rhrbqu

用ssega方法，结合GBM_ER文献里的mmune里来看

geneset = rio::import("mmc3.xlsx",skip = 2)
geneset = split(geneset$Metagene,geneset$`Cell type`)
lapply(geneset[1:3], head)
#得到每个免疫细胞有哪些基因型
## $`Activated B cell`
## [1] "ADAM28" "CD180"  "CD79B"  "BLK"    "CD19"   "MS4A1" 
## 
## $`Activated CD4 T cell`
## [1] "AIM2"  "BIRC3" "BRIP1" "CCL20" "CCL4"  "CCL5" 
## 
## $`Activated CD8 T cell`
## [1] "ADRM1"     "AHSA1"     "C1GALT1C1" "CCT6B"     "CD37"      "CD3D"


library(GSVA)
f = "ssgsea_result.Rdata"
if(!file.exists(f)){
  p = list()
for(i in 1:4){
  #i = 1
  re <- gsva(exps[[i]], geneset, method="ssgsea",
               mx.diff=FALSE, verbose=FALSE
)#得到的re，一行一个免疫细胞，一列一个样本
#做箱线图
  p[[i]] = draw_boxplot(re,survs[[i]]$Risk)+
    ggtitle(corhorts[[i]])+
    theme(plot.title = element_text(hjust = 0.5))+ 
    theme(plot.margin=unit(c(1,3,1,1),'lines')) #调整边距
}
  save(re,p,file = f)
}
load(f)

wrap_plots(p)+ plot_layout(guides = "collect") &
  theme(legend.position='right')

得到的re：

得到的可视化

3.泛癌免疫浸润的可视化，加强调色的泛癌亚型

https://www.yuque.com/xiaojiewanglezenmofenshen/bsgk2d/vhsx97ohpoxn4nwg

泛癌的亚型如何寻找，新版TCGA文献 ppt里的第23张

4.estimate算法

https://www.yuque.com/xiaojiewanglezenmofenshen/bsgk2d/yenun9

5.肿瘤突变负荷的计算

https://www.yuque.com/xiaojiewanglezenmofenshen/bsgk2d/gbyixw

maf文件，给每个病人计算他的TMB，为一个连续性数值，根据这个数值的大小把病人分成两个组，小于中位数的一个组，大于中位数的为另一个组

7.PCA图3D版

library(tinyarray)
?draw_pca
draw_pca(t(iris,1:4),iris$Species)
draw_pca(t(iris,1:4),iris$Species,style = "ggplot2")
draw_pca(t(iris,1:4),iris$Species,style = "3D")

8.突变频谱图

https://mp.weixin.qq.com/s/8ROtJxEXec9ukvty56QcbA

主要是数据下载难？

展示你想展示的那组基因的突变情况

options(stringsAsFactors = F) 
require(maftools) 
require(dplyr)
project='TCGA_KIRC'
laml = read.maf(maf = 'TCGA.KIRC.mutect.somatic.maf.gz')
laml 
#> An object of class  MAF 
#>                         ID summary   Mean Median
#>  1:             NCBI_Build      NA     NA     NA
#>  2:                 Center      NA     NA     NA
#>  3:                Samples     370     NA     NA
#>  4:                 nGenes   10243     NA     NA
#>  5:        Frame_Shift_Del    1792  4.843      4
#>  6:        Frame_Shift_Ins     463  1.251      1
#>  7:           In_Frame_Del     270  0.730      0
#>  8:           In_Frame_Ins      26  0.070      0
#>  9:      Missense_Mutation   17050 46.081     42
#> 10:      Nonsense_Mutation    1297  3.505      3
#> 11:       Nonstop_Mutation      33  0.089      0
#> 12:            Splice_Site     672  1.816      1
#> 13: Translation_Start_Site      34  0.092      0
#> 14:                  total   21637 58.478     53
length(unique(laml@data$Tumor_Sample_Barcode))
#> [1] 370
length(unique(laml@data$Hugo_Symbol))

需要的数据是maf，读取进来后，再用maftool可视化

plotmafSummary

9.WGCNA 加权基因共表达分析

weighted gene co-expression network analysis.寻找协同表的基因模块，探索基因网络与研究性状之间的联系。

每个表型相关模块里的那些基因

模块：具有高拓扑重叠相似性的基因合集。共表达模块是根据非相似性矩阵，利用聚类算法获得。基因与他所属的同一模块内的其他基因往往具有更高的共表达特性。

ME：代表模块的第一主分，即PCA1。用来描述模块在各样本中的表达模式。

MM：代表给定基因和模块ME之间的相关系数，描述基因属于一个模块的可靠性。该概念在模块划分时使用。

Hubhub基因代表强关联度的基因，往往有高的MM值

模块内连通性：某一基因的的模块内连通性同于该基因与模块内其他关联度之和，值越大说明这个基因在这个模块越处于核心位置

整体连通性等于给定基因和整个网络中其他基因关联度之和，描述该基因在网络中的地位

核心

关心模块与性状之间的关系→再找到hubgene，而非单个基因与性状的关系！

https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/ytog00

分析步骤

A.输入数据：表达矩阵

B.构建共表达网络

C.将表达模式相近的基因划分为一个模块

（模块划分➡合并相似模块)

D.模块与性状之间的关联分析，找到与目标性状相关性最高的模块，对相关性最高的模块的所有基因进行可视化展示（模块之间的关联分析）

从相关性最高的模块中筛选最重要的基因

E.模块中核心基因的鉴定（计算GS值（Gene significance）。找到GS值最高的基因，并对其进行qPCR分析。知道基因的功能，还有上游分析？（基因表观遗传学：甲基化水平与性状的关系）

F.得到结论

数据准备：基因的表达量、样本、每个样本的（某一个关心的）性状的表达量

单细胞之多样本整合

1.Harmoy整合多细胞数据

https://www.yuque.com/xiaojiewanglezenmofenshen/grgtd5/wfzdzp?singleDoc#

2.CCG+MNA（找到锚定细胞，去除批次效应

MNA(mutual nearest neighbor)算法以此找到两个数据集之间互相“距离”最近的细胞，Seurat将这些相互最近邻细胞称为“锚点细胞”（anchors）。

对应代码GSE117570

①对数据先进行标准化，并识别 variable feature

for (i in 1:length(scelist)) {
  scelist[[i]] <- NormalizeData(scelist[[i]], verbose = FALSE)#对数据进行标准化
  scelist[[i]] <- FindVariableFeatures(scelist[[i]],verbose = FALSE,nfeatures = 8000)#识别VariableFeatures
}

②找到Integration的基因，然后找到anchors细胞

features <- SelectIntegrationFeatures(object.list = scelist,nfeatures = 8000)#找到IntegrationFeatures
sce.anchors <- FindIntegrationAnchors(object.list = scelist,anchor.features = features)#利用IntegrationFeatures，并使用FindIntegrationAnchors函数识别锚点
sce.integrated <- IntegrateData(anchorset = sce.anchors)
#然后我们将这些锚点传递给IntegrateData函数，该函数返回一个Seurat对象。

#sce.integrated是三个样品经过批次效应矫正后合并的Seurat 对象，对这个对象进行分群分析

DefaultAssay(sce.integrated) <- "integrated"
sce.all = sce.integrated
#然后我们可以使用这个新的表达矩sce.all阵进行下游分析和可视化。

9.monocle

newCellDataSet() 创造对象

需要的数据：

1.表达矩阵exprs:行名为基因，列名为细胞编号。

2.细胞表型信息phenoData：第一列细胞编号，其他列是细胞的相关信息。

3.基因注释featureData：第一列是基因编号，其他列是基因的对应信息。

exp的要求它的列与phenoData的行数，它的行数与featureData的行数相互匹配。featureData至少要有一列”gene_short_name”, 就是基因的symbol

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)
# Set assay to RSEM estimated counts
assay(fluidigm) <- assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4] 
sample_ann <- as.data.frame(colData(fluidigm))

gene_ann <- data.frame(
 gene_short_name = row.names(ct), 
 row.names = row.names(ct)
)

pd <- new("AnnotatedDataFrame",
 data=sample_ann)
fd <- new("AnnotatedDataFrame",
 data=gene_ann)

sc_cds <- newCellDataSet(
 ct, 
 phenoData = pd,
 featureData =fd,
 expressionFamily = negbinomial.size(),
 lowerDetectionLimit=1)
sc_cds

创建后的对象含有：表达矩阵：rows as features (usually genes) and columns as cells

使用 featureData and phenoData 函数可以获取基因和样本信息

其中 expressionFamily指定表达矩阵的归一化形式

> HSMM<-monocle_cds
## 归一化 
> HSMM <- estimateSizeFactors(HSMM)
> HSMM <- estimateDispersions(HSMM)
#Filtering low-quality cells
> HSMM <- detectGenes(HSMM, min_expr = 3 )
> print(head(fData(HSMM)))#fData提取featureData
              gene_short_name num_cells_expressed
AL627309.1         AL627309.1                   9
AP006222.2         AP006222.2                   3
RP11-206L10.2   RP11-206L10.2                   5
RP11-206L10.9   RP11-206L10.9                   3
LINC00115           LINC00115                  18
NOC2L                   NOC2L                 254
> expressed_genes <- row.names(subset(fData(HSMM),
+                                     num_cells_expressed >= 10))
> print(head(pData(HSMM)))#pData提取pheoData
               orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters Size_Factor num_genes_expressed
AAACATACAACCAC     pbmc3k       2419          779  3.0177759               1               1          NA                 779
AAACATTGAGCTAC     pbmc3k       4903         1352  3.7935958               3               3          NA                1352
AAACATTGATCAGC     pbmc3k       3147         1129  0.8897363               1               1          NA                1129
AAACCGTGCTTCCG     pbmc3k       2639          960  1.7430845               2               2          NA                 960
AAACCGTGTATGCG     pbmc3k        980          521  1.2244898               6               6          NA                 521
AAACGCACTGGTAC     pbmc3k       2163          781  1.6643551               1               1          NA                 781

scater

创建SingleCellExperiment对象

需要的数据：exp(行是基因，列是细胞)

细胞的表达矩阵（注释信息可以是细胞名称，组织来源，收集样品的时间点，处理条件等等, 行是细胞，列是属性

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

r 语言

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

r 语言

#新TCGA

#文献复现里的几种算法

作者已关闭评论

0 条评论

热度