CNV拷贝数变异分析（GISTIC、maftools）

CNV拷贝数变异分析是什么？贴一段TCGA官网的介绍

"The copy number variation (CNV) pipeline uses Affymetrix SNP 6.0 array data to identify genomic regions that are repeated and infer the copy number of these repeats. This pipeline is built onto the existing TCGA level 2 data generated by Birdsuite and uses the DNAcopy R-package to perform a circular binary segmentation (CBS) analysis. CBS translates noisy intensity measurements into chromosomal regions of equal copy number. The final output files are segmented into genomic regions with the estimated copy number for each region. The GDC further transforms these copy number values into segment mean values, which are equal to log2(copy-number/ 2). Diploid regions will have a segment mean of zero, amplified regions will have positive values, and deletions will have negative values."

0. 目录

• 一. segment file数据下载和处理

• 1.1 从TCGA下载数据

• 1.2 数据处理

• 二. marker file数据下载和处理

• 2.1 从TCGA下载数据

• 2.2 提取freqcnv=FALSE数据，并且整理成标准格式

• 三. GenePattern GISTIC_2.0在线分析

• 四. maftools可视化GISTIC结果

1. segment file数据下载和处理

1.1 从TCGA下载数据

下载文件类型：
Copy Number Segment：A table that associates contiguous chromosomal segments with genomic coordinates, mean array intensity, and the number of probes that bind to each segment.

Masked Copy Number Segment：A table with the same information as the Copy Number Segment except that segments with probes known to contain germline mutations are removed

这里我用Masked Copy Number Segment做示范

rm(list = ls())
options(stringsAsFactors = F)
options(scipen = 200)

library(SummarizedExperiment)
library(TCGAbiolinks)

query <- GDCquery(project = "TCGA-BLCA",
                  data.category = "Copy Number Variation",
                  data.type = "Masked Copy Number Segment")
GDCdownload(query,method = "api")
BLCA_CNV_download <- GDCprepare(query = query, save = TRUE, save.filename = "BLCA_CNV_download.rda")

1.2数据处理

#读取rda文件
A=load("C:/Users/Meredith/Desktop/BLCA_CNV_download.rda")
tumorCNV <- eval(parse(text = A))

#改名
tumorCNV=tumorCNV[,2:7]
tumorCNV=tumorCNV[,c('Sample','Chromosome','Start','End','Num_Probes','Segment_Mean')]
write.table(tumorCNV,file = 'BLCA_CNV.txt',sep = '\t',quote = F,row.names = F)

BLCA_CNV.txt

#提取01A结尾的样本（这里我用了python，小伙伴们可以用R来做）
filename = 'BLCA_CNV.txt'
finalResultName = 'segment_file.txt'

read_file = open(filename)
out_file = open(finalResultName,"r+")
for line in read_file.readlines():
    data = line.split()
    x = data[0][13:16]
    if x == '01A':
        out_file.write(data[0])
        out_file.write('\t')
        out_file.write(data[1])
        out_file.write('\t')
        out_file.write(data[2])
        out_file.write('\t')
        out_file.write(data[3])
        out_file.write('\t')
        out_file.write(data[4])
        out_file.write('\t')
        out_file.write(data[5])
        out_file.write('\t')
        out_file.write('\n')

segment_file

2.marker file数据下载和处理

2.1从TCGA下载数据

TCGA现在的参考基因组版本是hg38，需要从官网下载marker file。下载地址：GDC Reference File Website，选择最新版本的“SNP6 GRCh38 Remapped Probeset File for Copy Number Variation Analysis”文件，并注意“If you are using Masked Copy Number Segment for GISTIC analysis, please only keep probesets with freqcnv =FALSE”

官网下载marker file

2.2提取freqcnv=FALSE数据，并且整理成标准格式

#这里我也是用的python（因为我电脑太菜用R要跑很久）
filename = "snp6.na35.remap.hg38.subset.txt"
finalResultName = "marker_file.txt"

read_file = open(filename)
out_file = open(finalResultName,"r+")
for line in read_file.readlines():
    data = line.split()
    if data[5]=='FALSE':
        out_file.write(data[0])
        out_file.write('\t')
        out_file.write(data[1])
        out_file.write('\t')
        out_file.write(data[2])
        out_file.write('\t')
        out_file.write('\n')

marker_file

3.GenePattern GISTIC_2.0在线分析

GenePattern

GISTIC_2.0
refgene file小伙伴们根据需要选择，这里我用的是TCGA下载的数据，所以选择hg38。将segment_file跟marker_file分别拖到seg file跟markers file区域。置信区间系统默认0.9，可以根据需要调整。点击RUN。
要等半个小时左右。
GISTIC_2.0输出结果

我们来看下其中两个文件

amp_qplot del_qplot

4.maftools可视化GISTIC结果

rm(list = ls())
options(stringsAsFactors = F)

BiocManager::install("PoisonAlien/maftools") #建议下载github最新版本的maftools包
library(maftools)

#读入GISTIC文件
laml.gistic = readGistic(gisticAllLesionsFile = 'C:/Users/Meredith/Desktop/maftools/all_lesions.conf_99.txt',
                         gisticAmpGenesFile = 'C:/Users/Meredith/Desktop/maftools/amp_genes.conf_99.txt', 
                         gisticDelGenesFile = 'C:/Users/Meredith/Desktop/maftools/del_genes.conf_99.txt', 
                         gisticScoresFile = 'C:/Users/Meredith/Desktop/maftools/scores.gistic', 
                         isTCGA = T)

readGistic里isTCGA=T作者给出的解释：please note that setting isTCGA=TRUE truncates sample names to first twelve characters - meaning primary samples and their metastatic counterparts (or from other sources) will all be collapsed into one single sample type. 所以在准备GISTIC输入文件时要将正常样品跟肿瘤样品分开。

##绘图
#genome plot
gisticChromPlot(gistic = laml.gistic, ref.build = 'hg38')

A

A. G-scores assigned by GISTIC for every cytoband plotted along the chromosome.

#Bubble plot
gisticBubblePlot(gistic = laml.gistic)

B

B. GISTIC results plotted as function of altered cytobands, mutated samples, and genes involved within the cytoband. Size of each bubble is according to -log10 transformed q values.

#oncoplot
gisticOncoPlot(gistic = laml.gistic, 
               sortByAnnotation = TRUE, top = 10)

C

C. Oncoplot displays most frequently altered (amplifications or deletions) copy number events ordered according to the frequency. Each columns represents a sample and each row represent a CNV segment.

#Visualizing CBS segments
plotCBSsegments(cbsFile = 'C:/Users/Meredith/Desktop/maftools/segment.file.txt', #这里的segment file的sample name只保留前12位
                tsb = 'ALL', #如果想指定某个样本，就输样本的名字，比如：tsb = 'TCGA-C4-A0F7'
                ref.build = 'hg38')

D

D. Plots segmented copy number data.