完成任意癌症的任意基因突变与否分组后的转录组测序的差异分析

[ixxmu]

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

[ixxmu]

完成任意癌症的任意基因突变与否分组后的转录组测序的差异分析 by 生信技能树

一个月前（ 2023-12-01 ）的学徒作业：任意癌症的任意基因突变与否分组后的转录组测序的差异分析，陆陆续续收到了一些反馈，有马拉松授课学员的也有学徒实习生的，发现虽然说给大家指明了数据分析结题思路，但大家仍然是千奇百怪的错误。总体上就5个步骤，大家可以错十几处：

• 从UCSC的XENA浏览器里面选择NSCLC的里面的LUAD数据集
• 然后下载LUAD的somatic的突变信息的maf文件
• 以及下载LUAD的表达量矩阵的counts文件
• 接着针对STK11基因把LUAD区分成为突变与否
• 最后根据STK11基因与否的分组信息进行差异分析，火山图和热图展示，富集分析等等

两种maf文件处理方式

如何找到somatic的突变信息的maf文件，仍然是从UCSC的XENA浏览器里面选择NSCLC的里面的LUAD数据集即可，这个是网页里面的鼠标点击操作。值得注意的是网页里面关于同一个癌症有两个跳转链接哦（其中一个带有GDC的前缀）：

• [TCGA Lung Adenocarcinoma (LUAD)](https://xenabrowser.net/datapages/?cohort=TCGA Lung Adenocarcinoma (LUAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (23 datasets)
• [GDC TCGA Lung Adenocarcinoma (LUAD)](https://xenabrowser.net/datapages/?cohort=GDC TCGA Lung Adenocarcinoma (LUAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (15 datasets)

首先可以看到不带GDC前缀的链接里面的突变主要是来源于mc3计划：

• https://tcga-xena-hub.s3.us-east-1.amazonaws.com/download/mc3_gene_level%2FLUAD_mc3_gene_level.txt.gz
• TCGA Unified Ensemble "MC3" gene-level mutation calls. 1: non-silent mutation 0: wt

MC3（Multi-Center Mutation Calling in Multiple Cancers）计划则是TCGA（The Cancer Genome Atlas）项目中的一个子项目，专注于对多个癌症类型进行突变信息的分析和整理。以下是MC3计划的主要特点和工作内容：

1. 突变信息整合： MC3计划旨在整合来自TCGA多个癌症类型的突变信息。这包括各种癌症样本中的单核苷酸变异（Single Nucleotide Variants，SNVs）、小型插入和缺失（Indels）等。
2. 标准化处理： 为了保证数据的一致性和可比性，MC3对从不同实验室和研究中心获得的数据进行了标准化处理。这涉及到对突变数据的质量控制、规范化和一致性检查。
3. 生物信息学方法： MC3采用了先进的生物信息学方法和工具，包括突变调用算法、数据处理流程等，以确保对基因组变异的准确和全面的分析。
4. 癌症类型涵盖： MC3覆盖了多个癌症类型，包括但不限于乳腺癌、肺癌、结肠癌等。这有助于研究者更全面地了解不同癌症类型的遗传变异。
5. 提供公共资源： MC3计划生成的数据被提供为公共资源，可以被科研机构、学者和生物信息学家免费使用。这种开放的数据共享有助于推动更广泛的癌症研究。
6. 突变谱分析： MC3对于突变的谱进行了分析，包括了基因的突变频率、突变类型的分布等，这对于了解癌症的遗传学机制和驱动突变具有重要意义。

有意思的是我读取这个MC3突变信息，发现本次作业涉及到的基因是没有的：

mut= data.table::fread('input/mc3_gene_level_LUAD_mc3_gene_level.txt.gz',data.table = F)
mut[1:4,1:4]
mut[mut$sample == 'STK11',]
mut[mut$sample == 'LKB1',]

> mut[1:4,1:4]
   sample TCGA-05-4244-01 TCGA-05-4249-01 TCGA-05-4250-01
1  UBE2Q2               0               0               0
2  CHMP1B               0               0               0
3 PSMA2P1               0               0               0
4  SHQ1P1               0               0               0

> mut$sample[which.max(rowSums(mut[,-1]))]
[1] "TP53"

大概率是因为这个MC3计划过于严格，过滤了绝大部分基因。不过突变病人数量最多的仍然是TP53基因，说明这个MC3信息本身是值得信赖的。

其次可以看到那个带GDC的链接进去就有4个不同的软件产出的somatic突变信息，如下所示 ：

4个不同的软件

我在生信技能树发布的很多 找somatic mutation教程大概率是都过时了，如下：

• TCGA计划的4个找somatic mutation的软件使用体验
• 肿瘤全外显子测序数据分析流程大放送
• 肿瘤空间异质性探究

毕竟是六年多过去了，然后在最新最全的mutect2教程，提到了其实大家不必在一棵树上吊死，而且如果是TCGA这样的公共数据集大概率是不需要自己从零开始处理fq文件后拿到somatic mutation信息。

直接下载即可，比如这个时候我们测试了muse这个软件的结果文件；

mut= data.table::fread('input/TCGA-LUAD.muse_snv.tsv.gz',data.table = F)
mut[1:4,1:4]
# somatic，4个软件，交集很少。。。
# 样品采集，送样，测序，到软件，参数，玄学。。。。
# maftools  
mutSid = mut$Sample_ID[mut$gene=='STK11']
head(mutSid)

> mut[1:4,1:4]
         Sample_ID     gene chrom    start
1 TCGA-35-4122-01A     AGRN  chr1  1046540
2 TCGA-35-4122-01A   PRDM16  chr1  3412324
3 TCGA-35-4122-01A PRAMEF11  chr1 12825031
4 TCGA-35-4122-01A PRAMEF20  chr1 13418296

> head(mutSid)
[1] "TCGA-44-7669-01A" "TCGA-44-7671-01A" 

可以看到的是已经是可以获取到那些有STK11基因突变的肺癌患者的id啦。

表达量矩阵的处理

仍然是从UCSC的XENA浏览器里面选择NSCLC的里面的LUAD数据集即可，值得注意的是选择带gdc前缀的哦，如下所示：

• log2(count+1)
• https://gdc-hub.s3.us-east-1.amazonaws.com/download/TCGA-LUAD.htseq_counts.tsv.gz

其实也是可以根据上面的网页链接的规律去获取所有的其它癌症的下载链接啦，然后就是读取TCGA-LUAD.htseq_counts.tsv.gz文件后的简单的处理，代码如下所示：

# 魔幻操作，一键清空
rm(list = ls()) 
options(stringsAsFactors = F)
library(data.table)
a1=fread( 'TCGA-LUAD.htseq_counts.tsv.gz' , 
          data.table = F) 
mat= a1[,2:ncol(a1)]  
mat=mat[1:(nrow(a1)-4),]
ensembl_matrix=ceiling(2^(mat)-1) #log2(x+1) transformed. 
rownames(ensembl_matrix) = gsub('[.][0-9]+','',a1$Ensembl_ID[1:(nrow(a1)-4)])
library(AnnoProbe) 
ids=annoGene(rownames(ensembl_matrix),'ENSEMBL','human') 
ids=ids[!duplicated(ids$SYMBOL),]
ids=ids[!duplicated(ids$ENSEMBL),]
symbol_matrix= ensembl_matrix[match(ids$ENSEMBL,
                                    rownames(ensembl_matrix)),] 
rownames(symbol_matrix) = ids$SYMBOL 
save(symbol_matrix,file = 'symbol_matrix.Rdata') 

分组后差异分析

需要根据突变信息对上面的表达量矩阵进行分组，所以是：

rm(list = ls()) 
library(data.table)
load(file = 'input/symbol_matrix.Rdata') 
group_list=ifelse(substring(colnames(symbol_matrix),14,15)=='11',
                  'control','case' )
table(group_list)
group_list = factor(group_list,levels = c('control','case' ))
symbol_matrix=symbol_matrix[,group_list=='case']

kp= !colnames(symbol_matrix) %in% unique(substring(mut$Sample_ID,1,15))
table(kp)
symbol_matrix=symbol_matrix[,kp]
group_list = ifelse(colnames(symbol_matrix) %in% mutSid,'case','control')
table(group_list) 
save(symbol_matrix,group_list,file = 'symbol_matrix.Rdata') 

有了表达量矩阵和分组信息，差异分析就是常规代码即可：

exprSet = symbol_matrix 
(colData <- data.frame(row.names=colnames(exprSet), 
                       group_list=group_list) )
dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = colData,
                              design = ~ group_list)
dds <- DESeq(dds) 
res <- results(dds, 
               contrast=c("group_list",
                          levels(group_list)[2],
                          levels(group_list)[1]))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG =as.data.frame(resOrdered)
DEG_deseq2 = na.omit(DEG)
   
save(DEG_deseq2, 
     file =  'DEG_deseq2.Rdata' ) 

大家赶快使用上面的代码去测试一下其它癌症吧，任意癌症的任意基因突变与否分组都可以。

其它看癌症突变全景图的方式

下面的网页链接里面的癌症的缩略词替换即可访问任意癌症的突变全景图：

• http://firebrowse.org/iCoMut/?cohort=LAML

但是突变本身有很多不同的概念，germline和somatic，snv和indel，而且snv还根据是否影响蛋白质细分很多不同肿瘤。

1. Germline 和 Somatic 变异：

• Germline 变异： 这些变异存在于生物个体的生殖细胞中（例如卵子和精子），并因此可以传递给后代。Germline 变异通常在个体的整个生命周期中存在，并影响其全部细胞。它们是个体遗传信息的一部分。
• Somatic 变异： 这些变异发生在个体的非生殖细胞（体细胞）中，只会影响个体的一部分细胞，而不会被传递给后代。Somatic 变异通常与癌症等疾病的发生有关，因为它们可能导致细胞生长和分裂的异常。

SNV 和 Indel：

• SNV（Single Nucleotide Variant）： 这是指单个核苷酸的变异，例如由于一个碱基被替换为另一个碱基而引起的变异。SNV 可以是点突变，包括单核苷酸替代、插入或缺失。
• Indel（Insertion/Deletion）： 这是指插入或删除一个或多个核苷酸的变异。Indel 可能导致基因框架的移位，影响蛋白质编码的读框，从而产生功能性影响。

蛋白质影响的 SNV 分类：

• Missense 变异： 单个氨基酸被替换为另一个氨基酸，这可能影响蛋白质功能。
• Nonsense 变异： 这种变异导致一个氨基酸被提前终止，导致蛋白质合成中止，通常会导致非功能性蛋白质的产生。
• Frameshift 变异： 插入或删除核苷酸导致读框发生改变，通常导致蛋白质合成的严重改变。
• Silent 变异： 核苷酸改变，但不会改变相应氨基酸的残基。

文末友情宣传

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

• 生物信息学马拉松授课（买一得五） ，你的生物信息学入门课，春节前最后一次
• 时隔5年，我们的生信技能树VIP学徒继续招生啦
  
• 144线程640Gb内存服务器共享一年仍然是仅需800