下載多個(gè)TCGA maf文件
每個(gè)腫瘤都有4個(gè)maf文件蟀俊,是用不同的軟件找出的突變信息钦铺,所以只要挑出一個(gè)就好。然后合并多個(gè)腫瘤的maf文件肢预,直接文件追加就行矛洞。
MutSigCV進(jìn)行突變負(fù)荷分析尋找Driver Gene
http://www.reibang.com/p/c86fadbff4c3
這個(gè)文章里講的很詳細(xì),生成的 sig_genes.txt 文件烫映,里邊有基因和p值沼本,代表基因?yàn)轵?qū)動(dòng)基因的可信度。
補(bǔ)充:TCGA里的maf文件對(duì)應(yīng)的參考基因組是GRCh38锭沟。所以不能使用MutSigCV官網(wǎng)的hg19抽兆,上面的文章里也寫(xiě)有方法,下面是實(shí)施的具體步驟冈钦。其中刪除換行符時(shí)候郊丛,用 sed -e ':a;N;s/\n//;ta' filename
沒(méi)有成功李请,似乎文件太大?所以使用下面方法厉熟。處理完后最好 wc 命令看下染色體長(zhǎng)度是不是官方標(biāo)明的長(zhǎng)度相等导盅,確保參考基因組構(gòu)建成功了。否則運(yùn)行程序的時(shí)候會(huì)有如下報(bào)錯(cuò)Error using MutSigCV>MutSig_preprocess (line 542) probable build mismatch between mutation_file and chr_files
誰(shuí)要是懶得構(gòu)建揍瑟,我發(fā)你郵箱白翻。
ls | while read file;do sed -i '1d' $file;done #刪除第一行
ls | while read file;do cat $file | tr "\n" " " >${file}new;done #新生成的文件里是換行符換成了空格
ls | while read file;do sed s/[[:space:]]//g $file > ${file}new;done #刪除空格
驅(qū)動(dòng)基因與患者信息對(duì)應(yīng)
根據(jù)p值篩選前5000個(gè)基因,然后轉(zhuǎn)置成0绢片,1(突變?yōu)?)滤馍。
library(data.table)
library(dplyr)
sig_genes=read.table(file = "output.sig_genes.txt",sep = "\t",header = T,stringsAsFactors=F)
maf=fread("my_input_data.maf",sep="\t",header=T,
fill=T,blank.lines.skip=T,stringsAsFactors=F)
genes=sig_genes$gene[1:5000]
samples=unique(maf$Tumor_Sample_Barcode)
result=matrix(data=NA,length(genes),length(samples))
n=1
for(s in samples){
g=which(maf$Tumor_Sample_Barcode==s) %>% maf$Hugo_Symbol[.]
l=genes %in% g
result[,n]=ifelse(l,1,0)
n=n+1
}
result.table=as.data.frame(result)
result.table.1=rbind(samples,result.table)
result.table.1=cbind(c("gene",genes),result.table.1)
write.table(result.table.1,file = "result_matrix.txt",
quote=F,sep="\t",row.names = F,col.names = F)