新版TCGA的突变SNP数据添加临床信息
文章目录
- 加载数据和R包
- 读取数据
今天给大家演示下如何用自己的数据完成maftools
的分析,主要是snp文件和临床信息的制作,其实很简单,但是网络上的教程都说的不清楚。
这次我们直接用之前TCGA-COAD和TCGA-READ合并后的数据演示,合并教程请看前一篇推文(公众号翻历史推文)
加载数据和R包
因为现在的TCGA数据库不能直接下载4种制作好的maf文件了,需要自己整理,如果你还不知道怎么整理,请看这篇内容(去翻推文)
load(file = "./TCGA-colrectum/colrectal_snp.rdata")library(maftools)
除此之外,我们还要给这个数据添加临床信息,也是只要先加载之前合并后的数据。
load(file = "./TCGA-colrectum/colrectal_clin.rdata")
maftools
包添加临床信息的方式非常简单,只要两个文件都有Tumor_Sample_Barcode
这一列且对得上就行。所以我们还要对snp文件和临床信息进行一些简单的处理。
- 对于两个文件中的
Tumor_Sample_Barcode
这一列,我们只要前12个字符即可 - 临床信息中有一些是Normal的样本,需要去除
- 之选择在snp文件中有的样本
# 只要前12个字符
colrec_snp$Tumor_Sample_Barcode <- substr(colrec_snp$Tumor_Sample_Barcode,1,12)
head(colrec_snp$Tumor_Sample_Barcode)
## [1] "TCGA-D5-6530" "TCGA-D5-6530" "TCGA-D5-6530" "TCGA-D5-6530" "TCGA-D5-6530"
## [6] "TCGA-D5-6530"
如果你是像我这样直接用的TCGAbiolinks
包下载的数据,那这个临床信息直接包含了sample_type
这一列,不需要自己根据样本名确定到底是normal还是tumor,十分方便。
index <- unique(colrec_snp$Tumor_Sample_Barcode)# 只要肿瘤样本
clin_snp <- clin[!clin$sample_type == "Solid Tissue Normal", ]# 只要snp文件中有的样本
clin_snp <- clin_snp[clin_snp$patient %in% index, ]# clin中没有Tumor_Sample_Barcode这一列,直接添加一列
clin_snp$Tumor_Sample_Barcode <- clin_snp$patient
这样两个需要的文件就制作好了。
colrec_snp[1:5,1:5]
## X1 Hugo_Symbol Entrez_Gene_Id Center NCBI_Build
## 1 1 AGRN 375790 BCM GRCh38
## 2 1 ACAP3 116983 BCM GRCh38
## 3 1 CALML6 163688 BCM GRCh38
## 4 1 PRKCZ 5590 BCM GRCh38
## 5 1 WRAP73 49856 BCM GRCh38
clin_snp[1:5,1:5]
## barcode patient
## TCGA-A6-5664-01A-21R-1839-07 TCGA-A6-5664-01A-21R-1839-07 TCGA-A6-5664
## TCGA-D5-6530-01A-11R-1723-07 TCGA-D5-6530-01A-11R-1723-07 TCGA-D5-6530
## TCGA-AA-3556-01A-01R-0821-07 TCGA-AA-3556-01A-01R-0821-07 TCGA-AA-3556
## TCGA-AA-3818-01A-01R-0905-07 TCGA-AA-3818-01A-01R-0905-07 TCGA-AA-3818
## TCGA-AA-3660-01A-01R-1723-07 TCGA-AA-3660-01A-01R-1723-07 TCGA-AA-3660
## sample shortLetterCode
## TCGA-A6-5664-01A-21R-1839-07 TCGA-A6-5664-01A TP
## TCGA-D5-6530-01A-11R-1723-07 TCGA-D5-6530-01A TP
## TCGA-AA-3556-01A-01R-0821-07 TCGA-AA-3556-01A TP
## TCGA-AA-3818-01A-01R-0905-07 TCGA-AA-3818-01A TP
## TCGA-AA-3660-01A-01R-1723-07 TCGA-AA-3660-01A TP
## definition
## TCGA-A6-5664-01A-21R-1839-07 Primary solid Tumor
## TCGA-D5-6530-01A-11R-1723-07 Primary solid Tumor
## TCGA-AA-3556-01A-01R-0821-07 Primary solid Tumor
## TCGA-AA-3818-01A-01R-0905-07 Primary solid Tumor
## TCGA-AA-3660-01A-01R-1723-07 Primary solid Tumor
读取数据
两个文件都没有问题,直接读取即可。
colrec.maf <- read.maf(colrec_snp,clinicalData = clin_snp)
## -Validating
## --Removed 41322 duplicated variants
## -Silent variants: 67679
## -Summarizing
## --Mutiple centers found
## BCM;WUGSC;BCM;WUGSC;BCM;BI--Possible FLAGS among top ten genes:
## TTN
## SYNE1
## MUC16
## -Processing clinical data
## --Annotation missing for below samples in MAF:
## TCGA-AA-3695
## TCGA-AA-3967
## TCGA-AF-2689
## TCGA-AF-3914
## TCGA-AG-3891
## TCGA-AG-3906
## TCGA-AZ-4681
## -Finished in 7.170s elapsed (6.480s cpu)
画图也是没有问题的,关于这个包的使用网络上教程非常多,也很全面,大家直接百度即可,这里就简单演示下。
plotmafSummary(colrec.maf)
最常见的额瀑布图:
oncoplot(maf = colrec.maf,clinicalFeatures = c("ajcc_pathologic_stage","vital_status"),top = 30,sortByAnnotation=T)
这个图其实就是Complexheatmap
画出来的,之前的推文里详细介绍了这个R包,感兴趣的可以看一看:
xxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxx
xxxxxxxxxx
colrec.titv = titv(maf = colrec.maf, plot = FALSE, useSyn = TRUE)
#plot titv summary
plotTiTv(res = colrec.titv)
比较有意思的棒棒糖图,这个图可以用trackviewer
包画的更好看,下次介绍。
画出来简单,但是解释过程挺复杂的,涉及到下游氨基酸的变化,需要了解碱基变化命名规则和氨基酸变化命名规则才能知道具体意思。
lollipopPlot(colrec.maf,gene = "TP53",AACol = "HGVSp_Short" # 需要注意这一列)
## 8 transcripts available. Use arguments refSeqID or proteinID to manually specify tx name.
## HGNC refseq.ID protein.ID aa.length
## 1: TP53 NM_000546 NP_000537 393
## 2: TP53 NM_001126112 NP_001119584 393
## 3: TP53 NM_001126118 NP_001119590 354
## 4: TP53 NM_001126115 NP_001119587 261
## 5: TP53 NM_001126113 NP_001119585 346
## 6: TP53 NM_001126117 NP_001119589 214
## 7: TP53 NM_001126114 NP_001119586 341
## 8: TP53 NM_001126116 NP_001119588 209
## Using longer transcript NM_000546 for now.
拷贝数变异肯定也是没有问题的,也是用的之前合并后的数据,然后经过gistic处理,就得到了我们需要的文件,关于gistic这个软件的使用,大家百度即可~
all.lesions <- "./TCGA-colrectum/TCGA_COREAD_results/all_lesions.conf_90.txt"
amp.genes <- "./TCGA-colrectum/TCGA_COREAD_results/amp_genes.conf_90.txt"
del.genes <- "./TCGA-colrectum/TCGA_COREAD_results/del_genes.conf_90.txt"
scores.gis <- "./TCGA-colrectum/TCGA_COREAD_results/scores.gistic"colrec.gistic = readGistic(gisticAllLesionsFile = all.lesions, gisticAmpGenesFile = amp.genes, gisticDelGenesFile = del.genes, gisticScoresFile = scores.gis, isTCGA = TRUE)
## -Processing Gistic files..
## --Processing amp_genes.conf_90.txt
## --Processing del_genes.conf_90.txt
## --Processing scores.gistic
## --Summarizing by samplescolrec.gistic
## An object of class GISTIC
## ID summary
## 1: Samples 611
## 2: nGenes 2791
## 3: cytoBands 88
## 4: Amp 97455
## 5: Del 313297
## 6: total 410752
结果也是毫无问题!
准备自己的数据,就是如此的简单。
新版TCGA的突变SNP数据添加临床信息相关推荐
- 新版TCGA不同癌种数据合并
很多文章对于TCGA中的一些癌症都是联合分析的,比如TCGA-COAD和TCGA-READ,首先是它们的疾病特点和治疗方式存在很多相似之处,同时这样做也可以增大样本量. 如果你是使用TCGAbioli ...
- R语言 关联TCGA数据库下载的RNA-SEQ数据和临床信息
刚开始学习TCGA数据处理和分析,记下来方便以后查看 setwd("E:/MyData/luadRNA-SEQ-20201028") #把工作目录定位到manifest文件所在的位 ...
- TCGAbiolinks整理表达数据和临床数据
新版TCGAbiolinks的整理表达数据和临床数据 没有废话直接干 ##加载包 rm(list = ls()) options(stringsAsFactors = F) gc() library( ...
- maftools|TCGA肿瘤突变数据的汇总,分析和可视化
之前介绍了使用maftools | 从头开始绘制发表级oncoplot(瀑布图) R-maftools包绘制组学突变结果(MAF)的oncoplot或者叫"瀑布图",以及一些细节的 ...
- 新版TCGA数据库学习:批量下载新版TCGA数据
众所周知,TCGA数据库改版了!!改的比之前更好用了! 对于常规转录组数据,主要是以下几点改变: 下载一次即可获得counts.TPM.FPKM三种类型的表达矩阵,再也不用单独下载了 自带gene s ...
- 最新版TCGA 矩阵整理,百分百复现成功
最近TCGA更新了,下载研究一下,我们从TCGA下载STAD的数据,选择其中的一个打开,发现了一个好消息那就是矩阵的整合难度降低了,而且提供TPM以及FPKM 还有校正的count 以及gene_na ...
- TCGA数据集介绍及数据下载
目录 一.TCGA数据集介绍 1.1 数据集介绍 1.2 File介绍 1.2.1 Data Category(数据类别) 1.2.2 Data Type(数据类型) 1.2.3 Experiment ...
- TCGA数据集介绍及数据下载指南(新手友好篇)
目录 一.TCGA数据集介绍 1.1 数据集介绍 1.2 File介绍 1.2.1 Data Category(数据类别) 1.2.2 Data Type(数据类型) 1.2.3 Experiment ...
- 区域卫生数据用于临床疗效分析的可用性研究
区域卫生数据用于临床疗效分析的可用性研究 叶琪1,赵亮1,阮彤1,冯东雷2,高炬3,刘珉3 1. 华东理工大学,上海 200237 2. 万达信息股份有限公司,上海 200040 3. 上海中医药大学 ...
- 新版TCGA表达矩阵1行代码提取2.0版
配合视频教程使用更佳:[1行代码提取6种TCGA表达矩阵和临床信息] https://www.bilibili.com/video/BV12R4y197Ne/?share_source=copy_we ...
最新文章
- 以太坊又一次大拥堵何去何从?深度对话美图以太坊DPoS算法实现团队
- Error: Cannot find module ‘webpack-cli/bin/config-yargs‘
- 开启Mysql慢查询来优化mysql
- Ubuntu 里的Spyder不能切换中文输入
- linux 文件服务,Linux操作系统之文件服务(ftp、nfs)
- Java笔记-Log4j在Spring Boot中的使用
- APUE 线程的分离状态
- ASP.NET 数据访问类
- 使用FragmentPagerAdapter和FragmentStatePagerAdapter时Fragment生命周期区别
- Python 列表(list)、字典(dict)、字符串(string)常用基本操作小结
- android NDK如何解决Please define the NDK_PROJECT_PATH variable to point to it
- 树莓派python 简介_树莓派与python语言概述
- 最新北京人才公寓申请流程,技术员的福利~
- 美国服务器怎么样 RAKsmart美国服务器适合做什么
- 机器学习算法——概率类模型评估指标4(校准可靠性曲线及预测概率直方图)
- 集成学习_GBDT_XGBoost
- 第二证券|沪指缩量跌0.25%,地产、医药板块强势,钠电池概念表现亮眼
- 库克逼腾讯分成30%遭拒,苹果APP Store或将微信下架!
- 【java】判断一个数是奇数还是偶数
- AIX小机不能启动故障-src:11002630代码故障的处理