科研信息学

TCGA 突变信息转置0,1 (机器学习输入)

2019-06-08  本文已影响41人  上校的猫

下载多个TCGA maf文件

每个肿瘤都有4个maf文件,是用不同的软件找出的突变信息,所以只要挑出一个就好。然后合并多个肿瘤的maf文件,直接文件追加就行。

MutSigCV进行突变负荷分析寻找Driver Gene

https://www.jianshu.com/p/c86fadbff4c3
这个文章里讲的很详细,生成的 sig_genes.txt 文件,里边有基因和p值,代表基因为驱动基因的可信度。
补充:TCGA里的maf文件对应的参考基因组是GRCh38。所以不能使用MutSigCV官网的hg19,上面的文章里也写有方法,下面是实施的具体步骤。其中删除换行符时候,用 sed -e ':a;N;s/\n//;ta' filename 没有成功,似乎文件太大?所以使用下面方法。处理完后最好 wc 命令看下染色体长度是不是官方标明的长度相等,确保参考基因组构建成功了。否则运行程序的时候会有如下报错Error using MutSigCV>MutSig_preprocess (line 542) probable build mismatch between mutation_file and chr_files
谁要是懒得构建,我发你邮箱。

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 #删除空格

驱动基因与患者信息对应

根据p值筛选前5000个基因,然后转置成0,1(突变为1)。

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)
上一篇下一篇

猜你喜欢

热点阅读