本文比较长,长到需要个目录 1.定义 2.TCGA拷贝数变异数据处理流程 2.1 数据下载 2.2 数据整理 3.GISTIC2.0 4.maftools可视化 找出感兴趣的区域里面的基因 5.基因层面的分析 5.1 文件和数值大小意义 5.2 差异CNV基因鉴定 5.3 棒棒糖图 6.与基因表达量联合分析 6.1 与表达量画箱线图 6.2 与表达量的相关性 7.生存分析
拷贝数变异(Copy number variation,CNV)是由基因组发生重排而导致的,一般指长度为1kb以上的基因组大片段的拷贝数增加或者减少,主要表现为亚显微水平的缺失和重复。CNV是基因组结构变异(Structuralvariation,SV)的重要组成部分。CNV位点的突变率远高于SNP(Single nucleotide polymorphism), 是人类疾病的重要致病因素之一。
–来自百度百科
TCGA 的拷贝数变异主要是Affymetrix SNP6.0 array,处理流程在: https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/CNV_Pipeline/
有三个workflow:ASCAT,ABSOLUTE和DNAcopy,其中ASCAT流程能够生成具有整数拷贝数值的等位基因特异性拷贝数段数据,以及派生的整数基因级拷贝数。参考基因组是hg19
ABSOLUTE 参考基因组是hg19,在GDC出版物页面有论文发表。
DNAcopy :进行循环二进制分割(CBS)将嘈杂的强度测量值转换为拷贝数相等的染色体区域。最终输出文件被分割成基因组区域,每个区域都有估计的拷贝数。GDC 进一步将这些拷贝数值转换为segment mean(段平均值),该值等于 log2(copy-number/ 2)。二倍体区域的段均值为零,扩增区域的区段均值为正值,缺失区的区段均值为负值。参考基因组是hg38,产生Copy Number Segment or Masked Copy Number Segment,其中Masked是指去掉了生殖相关SNV的数据。
我一般是喜欢用tcgabiolinks来下载,可以拿到最新的数据。
首先是下载Masked Copy Number Segment的数据,也就是上说的DNAcopy流程的结果。
rm(list = ls())
library(TCGAbiolinks)
f = "masked.Rdata"
if(!file.exists(f)){
query <- GDCquery(project = "TCGA-CHOL",
data.category = "Copy Number Variation",
data.type = "Masked Copy Number Segment")
GDCdownload(query)
seg <- GDCprepare(query)
colnames(seg)
seg = seg[,c(7,2:6)]
save(seg,file = f)
write.table(seg,file = "seg.txt",row.names = F,quote = F,sep = "\t")
}
load(f)
head(seg)
## # A tibble: 6 × 6
## Sample Chromosome Start End Num_Probes Segment_Mean
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 TCGA-W5-AA2X-11A-11D-A416-01 1 3.30e6 1.04e8 58516 -0.001
## 2 TCGA-W5-AA2X-11A-11D-A416-01 1 1.04e8 1.04e8 5 -0.972
## 3 TCGA-W5-AA2X-11A-11D-A416-01 1 1.04e8 1.79e8 27859 0.0016
## 4 TCGA-W5-AA2X-11A-11D-A416-01 1 1.79e8 1.79e8 2 -1.26
## 5 TCGA-W5-AA2X-11A-11D-A416-01 1 1.79e8 2.48e8 43396 0.0014
## 6 TCGA-W5-AA2X-11A-11D-A416-01 2 4.81e5 1.71e7 10742 0.008
length(unique(seg$Sample)) #人数
## [1] 85
#如果只需要tumor样本那么就:
#k = as.numeric(str_sub(seg$Sample,14,15))<10
#seg = seg[k,]
然后是下载markers file,在gdc网站的参考文件中能够找到它,下载后放在工作目录下。
https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files
此处有一句说明:If you are using Masked Copy Number Segment for GISTIC analysis, please only keep probesets with freqcnv = FALSE
marker = data.table::fread("snp6.na35.remap.hg38.subset.txt.gz",data.table = F)
head(marker)
## probeid chr pos strand type freqcnv par
## 1 SNP_A-1780270 7 78970267 + SNP TRUE FALSE
## 2 SNP_A-1780271 15 33103578 + SNP FALSE FALSE
## 3 SNP_A-1780272 1 189838554 - SNP TRUE FALSE
## 4 SNP_A-1780274 20 35320106 - SNP FALSE FALSE
## 5 SNP_A-1780277 12 75270266 + SNP FALSE FALSE
## 6 SNP_A-1780278 1 218717316 + SNP FALSE FALSE
可以看到数据里缺失有freqcnv 列,内容是TRUE和FALSE
table(marker$freqcnv)
##
## FALSE TRUE
## 1595984 241523
marker = marker[!marker$freqcnv,] #只保留FALSE
marker =marker[,1:3]
colnames(marker) = c("Marker Name","Chromosome","Marker Position")
head(marker)
## Marker Name Chromosome Marker Position
## 2 SNP_A-1780271 15 33103578
## 4 SNP_A-1780274 20 35320106
## 5 SNP_A-1780277 12 75270266
## 6 SNP_A-1780278 1 218717316
## 7 SNP_A-1780283 4 126709121
## 8 SNP_A-1780285 6 90209746
write.table(marker,file = "marker.txt",row.names = F,quote = F,sep = "\t")
colnames(marker)…那一句是把marker文件的列名修改为了GISTIC要求的名字。
至此,GISTIC2.0需要的输入文件就准备好了。
https://broadinstitute.github.io/gistic2/
全称是Genomic Identification of Significant Targets in Cancer(癌症中重要靶点的基因组鉴定),这是一个linux软件,GenePattern 网站上面有GISTIC模块,允许在线运行。
GISTIC可以识别基因组中在一组样本中被显著扩增或删除的区域,得到G-score和q值用于一组样本整体的可视化。
置信区间填0.99
然后拉到页面最后点run
正常的运行结果是有很多文件的,如果看到只有这三个文件那就是运行失败了。
出现这个结果的常见问题有:
1.markers文件要看清楚subset.txt.gz,不要下载错了 2.文件需要以tab为分隔符,所以上面的代码里sep = "\t"是不能少的。 3.markers文件的列名没改,所以上面的代码里colnames。。。那句代码是不能少的。
如果都排查完了还是三个文件,就试试不要上传markers。
打包下载并解压,放在工作目录下。
这个是区域层面的分析,单个基因层面的在下一部分。
library(maftools)
library(tidyverse)
laml.gistic = readGistic(gisticAllLesionsFile = "579436/all_lesions.conf_99.txt",
gisticAmpGenesFile = "579436/amp_genes.conf_99.txt",
gisticDelGenesFile = "579436/del_genes.conf_99.txt",
gisticScoresFile = "579436/scores.gistic",
isTCGA = TRUE)
## -Processing Gistic files..
## --Processing amp_genes.conf_99.txt
## --Processing del_genes.conf_99.txt
## --Processing scores.gistic
## --Summarizing by samples
head(laml.gistic@data)
## Hugo_Symbol Tumor_Sample_Barcode Variant_Classification Variant_Type
## 1: CCND1 TCGA-W5-AA30 Amp CNV
## 2: FGF4 TCGA-W5-AA30 Amp CNV
## 3: FGF19 TCGA-W5-AA30 Amp CNV
## 4: MYEOV TCGA-W5-AA30 Amp CNV
## 5: MRGPRD TCGA-W5-AA30 Amp CNV
## 6: MRGPRF TCGA-W5-AA30 Amp CNV
## Cytoband
## 1: AP_2:11q13.3
## 2: AP_2:11q13.3
## 3: AP_2:11q13.3
## 4: AP_2:11q13.3
## 5: AP_2:11q13.3
## 6: AP_2:11q13.3
ChromPlot,从左向右是全部的染色体。有颜色的部分突出了明显的扩增和删失区域,图上的文字例如”17q12”, “17” 指的是染色体编号,“q” 指的是长臂, “12” 是该臂上的具体区带或条带编号。
gisticChromPlot(gistic = laml.gistic, ref.build = 'hg38',markBands = "all")
head(laml.gistic@gis.scores)
## Variant_Classification Chromosome Start_Position End_Position fdr G_Score
## 1: Amp 1 3301765 3530249 0 0.003032
## 2: Amp 1 3558079 4117838 0 0.000000
## 3: Amp 1 4118243 4402844 0 0.003870
## 4: Amp 1 4403866 21803810 0 0.000000
## 5: Amp 1 21804530 22275547 0 0.004227
## 6: Amp 1 22276873 32428780 0 0.000000
## Avg_amplitude frequency
## 1: 0.406272 0.011765
## 2: 0.000000 0.000000
## 3: 0.258664 0.023529
## 4: 0.363623 0.011765
## 5: 0.142905 0.011765
## 6: 0.000000 0.000000
这个图的美化版本:https://blog.csdn.net/Ayue0616/article/details/127824394
Bubbleplot,横坐标是每个区域里发生拷贝数变异的样本数量,纵坐标是每个区域包含的基因数量,大小是按照每个区域的q值。
gisticBubblePlot(gistic = laml.gistic)
head(laml.gistic@cytoband.summary)
## Unique_Name nGenes nSamples Variant_Classification Cytoband
## 1: DP_1:1p36.21 110 31 Del 1p36.21
## 2: DP_13:9p21.3 2 23 Del 9p21.3
## 3: DP_5:4q34.3 69 23 Del 4q34.3
## 4: AP_2:11q13.3 13 10 Amp 11q13.3
## 5: DP_23:16q23.1 1 7 Del 16q23.1
## 6: DP_17:11q25 1635 14 Del 11q25
## Wide_Peak_Limits qvalues
## 1: chr1:8816334-14600926 4.9692e-26
## 2: chr9:21865499-21997723 2.7424e-07
## 3: chr4:176330306-190214555 4.6580e-07
## 4: chr11:68952543-69788268 2.4238e-04
## 5: chr16:78095161-79593873 2.0617e-03
## 6: chr11:1-135086622 3.5922e-03
OncoPlot,颜色太辣眼了必须得改改。。。
这个图一行是一个区域,一列是一个样本,展示每个区域在多少样本里发生了拷贝数变异。
mycolor = c("#f87669","#2fa1dd")
names(mycolor) = c("Amp","Del")
gisticOncoPlot(laml.gistic,top = 20,colors = mycolor)
tmp = as.data.frame(laml.gistic@data[,c(1,5)])
library(tidyverse)
tmp = separate(tmp,col = Cytoband,into = c("cnv","cyto"),sep = ":")
head(tmp)
## Hugo_Symbol cnv cyto
## 1 CCND1 AP_2 11q13.3
## 2 FGF4 AP_2 11q13.3
## 3 FGF19 AP_2 11q13.3
## 4 MYEOV AP_2 11q13.3
## 5 MRGPRD AP_2 11q13.3
## 6 MRGPRF AP_2 11q13.3
table(tmp$cnv)
##
## AP_1 AP_2 DP_1 DP_10 DP_11 DP_12 DP_13 DP_14 DP_15 DP_16 DP_17 DP_18 DP_19
## 12825 130 3410 5460 1944 20615 46 4452 18 5421 22890 12069 2227
## DP_2 DP_20 DP_21 DP_22 DP_23 DP_3 DP_4 DP_5 DP_6 DP_7 DP_8 DP_9
## 4760 28 28 6680 7 10136 17 1587 2168 150 26600 2222
k = tmp$Hugo_Symbol[tmp$cyto %in% c("1p36.21","11q25")] %>% unique()
length(k)
## [1] 1745
在GISTIC的输出结果里有4个genes.txt为结尾的文件
dir("579436/",pattern = "genes.txt$")
## [1] "all_data_by_genes.txt" "all_thresholded.by_genes.txt"
## [3] "broad_data_by_genes.txt" "focal_data_by_genes.txt"
其中:
all_data_by_genes.txt:包含基因在不同样本中具体的拷贝数数值。 all_thresholded.by_genes.txt:基因在不同样本中的拷贝数数值离散化后的结果,GISTIC2进一步将估计值阈值定为-2,-1,0,1,2,代表纯合子缺失、单拷贝缺失、二倍体正常拷贝、低水平拷贝数扩增或高拷贝数扩增。( homozygous deletion,single copy deletion,diploid normal copy,low-level copy number amplification,or high-level copy number amplification) broad_data_by_genes.txt:只考虑臂级(arm-level)事件的基因在不同样本中具体的拷贝数数值。 focal_data_by_genes.txt:只考虑局灶性(focal)事件的基因在不同样本中具体的拷贝数数值。
all_data = read.delim("579436/all_data_by_genes.txt",row.names = 1,check.names = F)
all_data[1:4,1:5]
## Gene ID Cytoband TCGA-W5-AA2X-11A-11D-A416-01
## ACAP3 116983 1p36.33 -0.003
## ACTRT2 140625 1p36.32 -0.003
## AGRN 375790 1p36.33 -0.003
## ANKRD65 441869 1p36.33 -0.003
## TCGA-W5-AA2X-10A-01D-A419-01 TCGA-W5-AA30-01A-31D-A416-01
## ACAP3 -0.001 -0.526
## ACTRT2 -0.001 -0.526
## AGRN -0.001 -0.526
## ANKRD65 -0.001 -0.526
range(all_data[,-(1:2)])
## [1] -1.293 3.657
all_thr = read.delim("579436/all_thresholded.by_genes.txt",row.names = 1,check.names = F)
all_thr[1:4,1:5]
## Locus ID Cytoband TCGA-W5-AA2X-11A-11D-A416-01
## ACAP3 116983 1p36.33 0
## ACTRT2 140625 1p36.32 0
## AGRN 375790 1p36.33 0
## ANKRD65 441869 1p36.33 0
## TCGA-W5-AA2X-10A-01D-A419-01 TCGA-W5-AA30-01A-31D-A416-01
## ACAP3 0 -1
## ACTRT2 0 -1
## AGRN 0 -1
## ANKRD65 0 -1
range(all_thr[,-(1:2)])
## [1] -2 2
tmp = all_thr[,-(1:2)]
table(tmp[tmp!=0])
##
## -2 -1 1 2
## 2021 218561 166880 8555
主要就看这两个文件啦。
对于这些数值大小的有好多种衍生的翻译
在timer2.0上面: “deep deletion”, “arm-level deletion”, “diploid/normal”, “arm-level gain”, and “high amplification”
在GSCA上面:Homozygous Deletion,Heterozygous Deletion,Heterozygous Amplification,Homozygous Deletion
Hete.(Heterozygous,杂合); Homo.(Homozygous,纯合)
翻了一些文献和github,可以根据是否有发生CNV将样本分成CNV(!=0)和Wild(==0),或者是分成扩增(>0)、缺失(<0)和正常(==0)三个组,看tumor和normal样本中每个基因的CNV是否有差别,检验方法主要是卡方检验。
Y叔团队开发的GeoTcgaData包有现成的函数differential_CNV,但我用的时候发现有个问题,源代码里面有一句gene <- gsub("\\..*", "", rownames(cnvData))
,代表把基因名中的点号后面的部分给去掉,这其实是给ensembel矩阵准备的。但如果是symbol为行名的矩阵,就要出错。
if(!require(GeoTcgaData))BiocManager::install("GeoTcgaData")
library(GeoTcgaData)
cnv = all_thr[,-(1:2)]
rownames(cnv) = 1:nrow(cnv)
cnv[cnv >= 1] = 1
cnv[cnv <= -1] = -1
library(tinyarray)
group = make_tcga_group(cnv)
diffcnv = differential_CNV(cnv,group)
对这个源代码进行修改就比较麻烦,不如把行名改成数字,然后再替换回symbol就是了。
diffcnv$gene = rownames(all_thr)
rownames(cnv) = rownames(all_thr)
head(diffcnv)
## gene P.Value adj.P.Val estimate
## 1 ACAP3 1.980127e-14 1.072074e-12 0
## 2 ACTRT2 1.980127e-14 1.072074e-12 0
## 3 AGRN 1.980127e-14 1.072074e-12 0
## 4 ANKRD65 1.980127e-14 1.072074e-12 0
## 5 ATAD3A 1.980127e-14 1.072074e-12 0
## 6 ATAD3B 1.980127e-14 1.072074e-12 0
table(diffcnv$P.Value<0.05)
##
## FALSE TRUE
## 1452 24536
table(diffcnv$adj.P.Val<0.05)
##
## FALSE TRUE
## 1452 24536
一组感兴趣的基因
my_gene = c("CXCL8", "FN1", "COL3A1", "ISG15", "COL1A2", "CXCL10", "ICAM1", "MMP9")
cnv = cnv[my_gene,group=="tumor"]
library(tidyverse)
cnv_df = as.data.frame(t(cnv)) %>% rownames_to_column("Sample")
cnv_long <- gather(cnv_df, key = "Gene", value = "CNV", 2:9)
# 统计每个基因的拷贝数变化样本比例
cnv_summary <- cnv_long %>%
filter(CNV != 0) %>% # 排除正常的样本
group_by(Gene, CNV) %>%
summarise(Freq = n()/ncol(cnv) *100) %>%
ungroup() %>%
mutate(CNV = ifelse(CNV>0,"Amp","Del"))
# 控制基因顺序
tmp = cnv_summary %>%
group_by(Gene) %>%
summarise(max = max(Freq)) %>%
arrange(desc(max))
cnv_summary$Gene = factor(cnv_summary$Gene,levels = tmp$Gene)
# 使用 ggplot2 绘制点图
library(ggplot2)
ggplot(cnv_summary, aes(x = Gene, y = Freq)) +
geom_segment( aes(x=Gene, xend=Gene, y=0, yend=Freq),color = "grey",linewidth = 3) +
geom_point(aes(color = CNV),size = 5) +
scale_color_manual(values = c("Del" = "#2fa1dd", "Amp" = "#f87669")) +
theme_minimal()
cnv_thr = as.matrix(all_thr[,colnames(cnv)])
cnvall = as.matrix(all_data[,colnames(cnv)])
library(tinyarray)
load("D:/TCGA_RNA_seq/tpm/chol_tpm.Rdata")
exp = trans_exp_new(chol_tpm)
exp = log2(exp+1)
colnames(exp) = str_sub(colnames(exp),1,16)
colnames(cnv_thr) = str_sub(colnames(cnv_thr),1,16)
colnames(cnvall) = str_sub(colnames(cnvall),1,16)
s = intersect(colnames(cnv_thr),colnames(exp))
exp = exp[my_gene,s]
cnv_thr = cnv_thr[,s]
cnvall = cnvall[,s]
g = "ISG15"
dat = data.frame(mrna = exp[g,],
cnv = as.factor(cnv_thr[g,]),
cnvscore = cnvall[g,])
head(dat)
## mrna cnv cnvscore
## TCGA-W5-AA30-01A 6.839125 -1 -0.526
## TCGA-W5-AA2W-01A 7.535962 -1 -0.889
## TCGA-ZH-A8Y4-01A 5.859495 -1 -0.802
## TCGA-W5-AA2Z-01A 8.296128 0 0.027
## TCGA-W6-AA0S-01A 7.077692 -1 -0.765
## TCGA-W5-AA33-01A 6.458971 -1 -0.409
table(dat$cnv)
##
## -2 -1 0 1
## 1 27 6 1
library(ggplot2)
ggplot(dat,aes(cnv,mrna,fill = cnv))+
geom_boxplot()+
scale_fill_brewer(palette = "Set1")+
theme_minimal()
library(ggstatsplot)
ggbetweenstats(dat,x = "cnv",y = "mrna",drop = F)
还是上面的数据
corscatterplot(dat,x = "cnvscore",y = "mrna")
surv = read.delim("D:/TCGA_RNA_seq/survival/TCGA-CHOL.survival.tsv",row.names = 1)
head(surv)
## OS X_PATIENT OS.time
## TCGA-W5-AA31-01A 0 TCGA-W5-AA31 10
## TCGA-W5-AA31-11A 0 TCGA-W5-AA31 10
## TCGA-WD-A7RX-01A 1 TCGA-WD-A7RX 21
## TCGA-YR-A95A-01A 1 TCGA-YR-A95A 26
## TCGA-4G-AAZR-11A 1 TCGA-4G-AAZR 28
## TCGA-4G-AAZR-01A 1 TCGA-4G-AAZR 28
colnames(surv)[c(1,3)] = c("event","time")
s = intersect(colnames(cnv_thr),rownames(surv))
surv = surv[s,]
cnv_thr = cnv_thr[,s]
g = "ISG15"
cnv_group = ifelse(cnv_thr[g,]==0,"Wild","CNV")
cnv_group = factor(cnv_group,levels = c("Wild","CNV"))
draw_KM(surv,cnv_group)
扫码关注腾讯云开发者
领取腾讯云代金券
Copyright © 2013 - 2025 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号:粤B2-20090059 深公网安备号 44030502008569
腾讯云计算(北京)有限责任公司 京ICP证150476号 | 京ICP备11018762号 | 京公网安备号11010802020287
Copyright © 2013 - 2025 Tencent Cloud.
All Rights Reserved. 腾讯云 版权所有