zhuang_gj
3月-09-2021
从TCGA下载Masked Copy Number Segment 数据-----解压后进行处理
rm(list = ls())
## 将整理好的Masked Copy Number Segment重新命名为“rawdata”
## 查看rawdata下面的文件
length(dir("rawdata") )
## [1] 1025
## 创建一个文件用来放rawdata里面的txt文件
dir.create("rawdata_all")
## 用lapply将rawdata里面txt文件整理到一个文件夹里
raw_txt <- dir("rawdata/")
lapply(raw_txt, function(x){
mydir <- paste0("./rawdata/",x)
file <- list.files(mydir,pattern = "nocnv_grch38")
myfile <- paste0("./rawdata/",x,"/",file)
file.copy(myfile,"rawdata_all")
})
## [[1]]
## [1] FALSE
##
## [[2]]
## [1] FALSE
library(data.table)
cnv_file <- dir("rawdata_all")
cnv_file[1]
## [1] "AMAZE_p_TCGASNP_b86_87_88_N_GenomeWideSNP_6_A01_735406.nocnv_grch38.seg.v2.txt"
## 文件所在位置
file <- paste0("./rawdata_all/",cnv_file)
head(file )
## [1] "./rawdata_all/AMAZE_p_TCGASNP_b86_87_88_N_GenomeWideSNP_6_A01_735406.nocnv_grch38.seg.v2.txt"
## [2] "./rawdata_all/AMAZE_p_TCGASNP_b86_87_88_N_GenomeWideSNP_6_A02_735476.nocnv_grch38.seg.v2.txt"
## [3] "./rawdata_all/AMAZE_p_TCGASNP_b86_87_88_N_GenomeWideSNP_6_A03_735420.nocnv_grch38.seg.v2.txt"
## [4] "./rawdata_all/AMAZE_p_TCGASNP_b86_87_88_N_GenomeWideSNP_6_A04_735430.nocnv_grch38.seg.v2.txt"
## [5] "./rawdata_all/AMAZE_p_TCGASNP_b86_87_88_N_GenomeWideSNP_6_A05_735540.nocnv_grch38.seg.v2.txt"
## [6] "./rawdata_all/AMAZE_p_TCGASNP_b86_87_88_N_GenomeWideSNP_6_A06_735442.nocnv_grch38.seg.v2.txt"
## 批量读取txt文件
file_list <- list()
fred_cnv <- function(x){
file_list[[x]] <- data.table::fread(file = x,data.table = F)
file_list[[x]]
}
cnv_df <- lapply(file, fred_cnv)
cnv_df <- do.call(rbind,cnv_df)
head(cnv_df)
## GDC_Aliquot Chromosome Start End
## 1 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 1 3301765 247650984
## 2 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 480597 236626512
## 3 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 236626820 236627088
## 4 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 236631315 237489539
## 5 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 237489625 237489879
## 6 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 237492059 239533237
## Num_Probes Segment_Mean
## 1 129636 0.0174
## 2 129172 0.0159
## 3 2 -1.5617
## 4 639 0.0003
## 5 3 -1.4417
## 6 1100 0.0153
## 读取metadata里面的注释信息
metadata <- jsonlite::fromJSON("metadata.cart.2021-02-15.json")
library(dplyr)
metadata_id <- metadata %>%
dplyr::select(c(file_name,associated_entities))
meta_df <- do.call(rbind,metadata$associated_entities)
head(metadata_id,2)
## file_name
## 1 TRIAL_p_TCGA_191_193_SNP_N_GenomeWideSNP_6_F12_932404.nocnv_grch38.seg.v2.txt
## 2 COAPT_p_TCGASNP_197_201_204_N_GenomeWideSNP_6_D01_1051508.nocnv_grch38.seg.v2.txt
## associated_entities
## 1 TCGA-EL-A3D5-10A-01D-A201-01, aliquot, e9ccf355-e68a-489c-a914-31b766bd0aef, a1e0d8be-7916-4bb2-9c3f-1e74b51c6ed3
## 2 TCGA-E3-A3DY-10A-01D-A20A-01, aliquot, dffef9e8-5b69-43d9-a95b-0744daa38d1d, b75decce-46f0-4689-be1b-08a9943775eb
## 匹配样本
cnv_df$Sample <- meta_df$entity_submitter_id[match(cnv_df$GDC_Aliquot,meta_df$entity_id)]
length(unique(cnv_df$GDC_Aliquot))
## [1] 1025
length(unique(cnv_df$Sample))
## [1] 1025
head(cnv_df)
## GDC_Aliquot Chromosome Start End
## 1 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 1 3301765 247650984
## 2 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 480597 236626512
## 3 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 236626820 236627088
## 4 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 236631315 237489539
## 5 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 237489625 237489879
## 6 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 237492059 239533237
## Num_Probes Segment_Mean Sample
## 1 129636 0.0174 TCGA-BJ-A0Z9-01A-11D-A10T-01
## 2 129172 0.0159 TCGA-BJ-A0Z9-01A-11D-A10T-01
## 3 2 -1.5617 TCGA-BJ-A0Z9-01A-11D-A10T-01
## 4 639 0.0003 TCGA-BJ-A0Z9-01A-11D-A10T-01
## 5 3 -1.4417 TCGA-BJ-A0Z9-01A-11D-A10T-01
## 6 1100 0.0153 TCGA-BJ-A0Z9-01A-11D-A10T-01
## 改名
cnv_df$Sample <- substring(cnv_df$Sample,1,16)
head(cnv_df)
## GDC_Aliquot Chromosome Start End
## 1 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 1 3301765 247650984
## 2 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 480597 236626512
## 3 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 236626820 236627088
## 4 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 236631315 237489539
## 5 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 237489625 237489879
## 6 3c7ea150-a119-4cd1-a7ef-1543ef87eebf 2 237492059 239533237
## Num_Probes Segment_Mean Sample
## 1 129636 0.0174 TCGA-BJ-A0Z9-01A
## 2 129172 0.0159 TCGA-BJ-A0Z9-01A
## 3 2 -1.5617 TCGA-BJ-A0Z9-01A
## 4 639 0.0003 TCGA-BJ-A0Z9-01A
## 5 3 -1.4417 TCGA-BJ-A0Z9-01A
## 6 1100 0.0153 TCGA-BJ-A0Z9-01A
## Gistic只需要maskedCNS的六列
cnv_df <- cnv_df[,c('Sample','Chromosome','Start','End','Num_Probes','Segment_Mean')]
head(cnv_df)
## Sample Chromosome Start End Num_Probes Segment_Mean
## 1 TCGA-BJ-A0Z9-01A 1 3301765 247650984 129636 0.0174
## 2 TCGA-BJ-A0Z9-01A 2 480597 236626512 129172 0.0159
## 3 TCGA-BJ-A0Z9-01A 2 236626820 236627088 2 -1.5617
## 4 TCGA-BJ-A0Z9-01A 2 236631315 237489539 639 0.0003
## 5 TCGA-BJ-A0Z9-01A 2 237489625 237489879 3 -1.4417
## 6 TCGA-BJ-A0Z9-01A 2 237492059 239533237 1100 0.0153
# 只挑选肿瘤样本
dim(cnv_df)
## [1] 61626 6
cnv_df <- cnv_df[grep("01A$",cnv_df$Sample),]
dim(cnv_df)
## [1] 28628 6
head(cnv_df)
## Sample Chromosome Start End Num_Probes Segment_Mean
## 1 TCGA-BJ-A0Z9-01A 1 3301765 247650984 129636 0.0174
## 2 TCGA-BJ-A0Z9-01A 2 480597 236626512 129172 0.0159
## 3 TCGA-BJ-A0Z9-01A 2 236626820 236627088 2 -1.5617
## 4 TCGA-BJ-A0Z9-01A 2 236631315 237489539 639 0.0003
## 5 TCGA-BJ-A0Z9-01A 2 237489625 237489879 3 -1.4417
## 6 TCGA-BJ-A0Z9-01A 2 237492059 239533237 1100 0.0153
## 将整理好的数据写出去
write.table(cnv_df,"MaskedCopyNumberSegment(Tumor).txt",sep="\t",
quote = F,col.names = F,row.names = F)
marker file数据下载和处理
GDC Reference Files---选择最新版本的“SNP6 GRCh38 Remapped Probeset File for Copy Number Variation Analysis”文件snp6.na35.remap.hg38.subset.txt.gz
注意：If you are using Masked Copy Number Segment for GISTIC analysis, please only keep probesets with freqcnv = FALSE
hg38_marker_file <- read.delim("snp6.na35.remap.hg38.subset.txt.gz")
head(hg38_marker_file )
## 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
# 注意“If you are using Masked Copy Number Segment for GISTIC analysis, please only keep probesets with freqcnv =FALSE”
hg_marker_file <- hg38_marker_file[hg38_marker_file$freqcnv=="FALSE",]
hg_marker_file <- hg_marker_file [,c(1,2,3)]
head(hg_marker_file )
## probeid chr pos
## 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(hg_marker_file,file = "hg_marker_file.txt",sep = "\t",col.names = T,row.names = F)
有兴趣的可以试试 TCGAbiolinks包看看结果如何。
下回讲讲如何用整理好的文件使用Gistic 2.0在线分析的得到Amplification GISTIC plot 以及Deletion GISTIC plot ~~~~~~~~~~~~~ maftools可视化相关结果 ~~~~~~~~~~~~~ 挑选拷贝数差异区域的基因。
坚持写代码的快一年半了，零零碎碎的代码产生的数据，占了我电脑大部分的空间，是时候整理下代码，让它变得更方便查询，毕竟每次写代码总会有改变，而有些是不需要变的。要感谢导师罗老师给予物质和精神上的支持，让我学习生物信息学，不给我设限，自由的探索数据。
