library(data.table)
library(foreach)
library(doMC)
library(ABSOLUTE)
setwd("/home/train/data/ABSOLUTE_test/Tumor_purity")
DoAbsolute <- function(scan,input) {
# registerDoSEQ() #多核心计算时候,需要注释掉
library(ABSOLUTE)
plate.name <- "test" #可以根据需要修改
genome <- "hg18"
platform <- c("SNP_6.0", "Illumina_WES", "SNP_250K_STY")
sample.name <- scan
primary.disease = "NA"
sigma.p <- 0
max.sigma.h <- 0.02
min.ploidy <- 0.95
max.ploidy <- 10
max.as.seg.count <- 1500
max.non.clonal <- 0
max.neg.genome <- 0
copy_num_type <- "total"
seg.dat.fn <- input
results.dir <- file.path(".", "output", scan, "absolute")
print(paste("Starting scan", scan, "at", results.dir))
log.dir <- file.path(".", "output", "abs_logs")
if (!file.exists(log.dir)) {
dir.create(log.dir, recursive=TRUE)
}
dz <- lapply(results.dir,function(x){
if(!file.exists(x)){
dir.create(x, recursive=TRUE)
}
})
rm(dz)
sink(file=file.path(log.dir, paste(scan, ".abs.out.txt", sep="")))
RunAbsolute(seg.dat.fn, sigma.p, max.sigma.h, min.ploidy, max.ploidy, primary.disease,
platform, sample.name, results.dir, max.as.seg.count, max.non.clonal,
max.neg.genome, copy_num_type, verbose=TRUE)
sink()
}
registerDoMC(20) #对每组数据调用线程池
arrays.txt <- "SNP6_solid_tumor.seg.txt"
input_data <- "test.txt"
scans <- unique(as.data.frame(fread(arrays.txt))$Sample)
foreach (scan=scans, .combine=c) %dopar% {
DoAbsolute(scan,input=input_data)
}
obj.name <- "test_summary"
results.dir <- file.path(".", "output", "abs_summary")
absolute.files <- file.path(".", "output",
scans, "absolute",
paste(scans, ".ABSOLUTE.RData", sep=""))
CreateReviewObject(obj.name, absolute.files, results.dir, "total", verbose=TRUE)
calls.path = file.path("output", "abs_summary", "test_summary.PP-calls_tab.txt")
modes.path = file.path("output", "abs_summary", "test_summary.PP-modes.data.RData")
output.path = file.path("output", "abs_extract")
ExtractReviewedResults(calls.path, "test", modes.path, output.path, "absolute", "total")
# doabsolute包更好使用,需要根据输入数据调整格式暂时分步处理,今后有时间写成pipline可以无缝对接doabsolute进行批量处理
mkdir T1_5
mkdir total
mv T1_5.cr.igv.seg T1_5
mv T1_5_absolute.maf T1_5
cd T1_5/
mv T1_5.cr.igv.seg tumor.seg
mv ../T1_3/step1_Xto23.sh .
bash step1_Xto23.sh
cp ../T1_3/step2_MafFilebalance.pl .
perl step2_MafFilebalance.pl T1_5_absolute.maf
cp ../T1_3/step3_SegFilebalance.pl .
perl step3_SegFilebalance.pl new_tumor.seg
cp ../T1_3/step4_rebalance_seg_and_maf.R .
Rscript step4_rebalance_seg_and_maf.R
cp ../T1_3/step5_get_seperate_file.R .
Rscript step5_get_seperate_file.R
cp ../T1_3/step6_prepare_doabsolute.pl .
perl step6_prepare_doabsolute.pl ./seperate_file/
cp ../T1_3/step6.5_prepare_doabsolute.sh .
#修改step6.5中的文件名为T1_5
bash step6.5_prepare_doabsolute.sh
cd ..
mkdir total
cd total/
mkdir seperate_file
cp ../T1_3/seperate_file/*.maf seperate_file/
cp ../T1_3/seperate_file/new_T1_3.seg seperate_file/
cp ../T1_4/seperate_file/*.maf seperate_file/
cp ../T1_4/seperate_file/new_T1_4.seg seperate_file/
cp ../T1_5/seperate_file/*.maf seperate_file/
cp ../T1_5/seperate_file/new_T1_5.seg seperate_file/
cp ../T1_3/new_final.seg .
cat ../T1_4/new_final.seg | awk 'NR>1' >> new_final.seg
cat ../T1_5/new_final.seg | awk 'NR>1' >> new_final.seg
cp ../T1_3/step7_doabsolute.R .
#注意修改路径
Rscript step7_doabsolute.R