CONIPHER: https:\/\/github.com\/McGranahanLab\/CONIPHER<\/a> conipher:<\/p>\n alpaca: \u9700gurobi\u8bc1\u4e66\uff0c\u672c\u5730\u5b89\u88c5<\/p>\n ASCAT: \u8fd0\u884calpaca\u9700\u8981refphase\uff0crefphase\u9700\u8981ASCAT<\/p>\n refphase:<\/p>\n \u5bf9\u591a\u7076\u6837\u672c\u8fd0\u884cMutect2\uff0cvarscan2\u540e\uff0c\u7b5b\u9009PASS\uff0c\u4f7f\u7528Annovar\u3001vcf2maf\u6ce8\u91ca\uff0c\u4e4b\u540e\u5206\u522b\u5408\u5e76\u6837\u672c\uff1a<\/p>\n \u5728R\u4e2d\u7b5b\u9009\u8fc7\u6ee4\uff0c\u53ea\u7559\u4e0b\u4e24\u4e2a\u8f6f\u4ef6\u5171\u540c\u9274\u5b9a\u5230\u7684\u3001gnomAD EAS (east asian population) AF<= 0.05\u30011000GenomesProject EAS <= 0.05, Tumor and normal total coverage >=20 & <=1000\uff0cTumor ALT read count >= 3\uff0cNormal ALT read count <= 1\uff0cTumor VAF >=0.05\uff0cNormal VAF <=0.05 \u7684\u7a81\u53d8<\/p>\n \u76f4\u63a5\u5220\u9664-BNPT\u548c-FA\uff0c\u83b7\u53d6\u5355\u72ec\u4ee5-M\u548c-E\u533a\u5206\u7684Tumor_Sample_Barcode\uff0c\u7559\u4e0b\u5982FETB01-M\u8fd9\u6837\u7684\u7f16\u53f7\uff0c\u8fd9\u6837\u53ef\u4ee5\u63d0\u53d6\u914d\u5bf9\u4e0a\u76ae\u3001\u57fa\u8d28\u6837\u672c\u4e2d\u7684PASS\u4f4d\u70b9<\/p>\n \u4f7f\u7528maf2vcf\u628a\u63d0\u53d6\u5b8c\u7684\u4f4d\u70b9\u8f6c\u6362\u4e3avcf\u6587\u4ef6\uff0c\u5176\u4e2d--per-tn-vcfs\u4f1a\u6839\u636eTumor_Sample_Barcode\uff0cvcf\u62c6\u5206\u4e3a\u6bcf\u4e2a\u914d\u5bf9\u6837\u672c\u4e00\u4e2a\uff08\u5982FETB01-M\u7684\u7a81\u53d8\u4f4d\u70b9\u4f1a\u5408\u5728\u4e00\u8d77\u6210\u4e00\u4e2avcf\uff09,\u751f\u6210\u7684vcf\u6587\u4ef6\u7684\u540d\u5b57\u4f1a\u662f\u201cFETB01-M_vs_FETB01-B.vcf\u201d\u7684\u683c\u5f0f<\/p>\n \u4e3acall Somatic Mutation \u4e2d\u7684 sample_pair\u6dfb\u52a0\u4e00\u5217(\u6700\u53f3\u8fb9)\uff0c\u6309\u7ec4\u7ec7\u7c7b\u578b\u533a\u5206\uff08FETB01-M\uff09<\/p>\n \u4f7f\u7528mutect2 force call \u6a21\u5f0f\uff0c\u4ee5\u4e0a\u4e00\u6b65\u751f\u6210\u7684vcf\u6587\u4ef6\u4e3atarget region\uff0ccall Somatic mutation\uff0c\u8fd9\u4e00\u6b65call\u5b8c\u7a81\u53d8\u540e\u4e00\u6837\u7528Annovar\u6ce8\u91ca\u3001\u5e76\u901a\u8fc7vcf2maf\u518d\u6b21\u8f6c\u6362\u4e3amaf\u6587\u4ef6\uff0c<\/p>\n \u81f3\u6b64\uff0c\u5236\u4f5c\u597d\u4e86\u7ecf\u8fc7\u8d28\u91cf\u8fc7\u6ee4\u3001\u6837\u672c\u914d\u5bf9\u7684\u7a81\u53d8\u7684maf\u6587\u4ef6 \u7531\u4e8econipher\u548crephase\u4e2d\uff0c\u6837\u672c\u540d\u5b57\u683c\u5f0f\u4e0d\u540c\uff08CONIPHER\u4e0d\u5141\u8bb8\u201c-\u201d\uff0c\u53ea\u80fd\u7528\u201c_\u201d\uff0c\u56e0\u6b64\u9700\u8981\u5728R\u4e2d\u8fdb\u884c\u4fee\u6539\uff09<\/p>\n \u8fd0\u884cALPACA<\/p>\n 0. \u5b98\u7f51\uff1a CONIPHER: https:\/\/github.com\/McGranahanLab\/CONIP…<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-619","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/619","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/comments?post=619"}],"version-history":[{"count":2,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/619\/revisions"}],"predecessor-version":[{"id":621,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/619\/revisions\/621"}],"wp:attachment":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/media?parent=619"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/categories?post=619"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/tags?post=619"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}
\nALPACA: https:\/\/github.com\/McGranahanLab\/ALPACA-model<\/a>
\nASCAT: https:\/\/github.com\/VanLoo-lab\/ascat<\/a>
\nrefphase: https:\/\/bitbucket.org\/schwarzlab\/refphase\/src\/master\/<\/a><\/p>\n1. Dependency\u8f6f\u4ef6\u5b89\u88c5<\/h2>\n
conda create -n conipher -c conda-forge -c bioconda conipher<\/code><\/pre>\ngit clone https:\/\/github.com\/McGranahanLab\/ALPACA-model.git\ncd ALPACA-model\nconda env create --name alpaca --file environment.yml\nconda run -n alpaca pip install dist\/*.whl<\/code><\/pre>\n# R\nBiocManager::install(c('GenomicRanges','IRanges')\ndevtools::install_bitbucket("schwarzlab\/ascat_v3_fork\/ASCAT") # \u9700\u8981\u5b89\u88c5\u4ed6\u4eec\u4fee\u6539\u8fc7\u7684ASCAT V3\u624d\u884c<\/code><\/pre>\nBiocManager::install("GenomicRanges")\ndevtools::install_bitbucket("schwarzlab\/refphase")<\/code><\/pre>\n2. CONIPHER<\/h2>\n
2.1 \u51c6\u5907\u8f93\u5165\u6587\u4ef6<\/h3>\n
# mutect2\ncat *.hg38_multianno.maf | awk 'NR<=3 || (!\/^#\/ && $1!="Hugo_Symbol")' > 6patients_mutect2_pass.Somatic.hc.hg38_multianno.maf\n\n# varscan2\ncat *.hc.hg38_multianno.maf | awk 'NR<=3 || (!\/^#\/ && $1!="Hugo_Symbol")' > 6patients_varscan2_pass.Somatic.hc.hg38_multianno.maf<\/code><\/pre>\nmutect2_maf_df <- read.table("\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/Mutation\/mutect2_TNpaired\/6patients_mutect2_pass.Somatic.hc.hg38_multianno.maf", header = T, sep = '\\t', quote = "")\nvarscan2_maf_df <- read.table("\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/Mutation\/varscan2\/6patients_varscan2_pass.Somatic.hc.hg38_multianno.maf", header = T, sep = '\\t', quote = "")\n\ntable(mutect2_maf_df$Tumor_Sample_Barcode)\ntable(varscan2_maf_df$Tumor_Sample_Barcode)\n\nfilter_maf_df <- function(maf_df){\n # reformat gnomAD EAS (east asian population) AF and only keep AF <= 0.05\n maf_df$gnomADe_EAS_AF[is.na(maf_df$gnomADe_EAS_AF)] <- 0\n maf_df <- subset(maf_df,gnomADe_EAS_AF <= 0.05)\n\n # reformat 1000GenomesProject EAS (east asian population) AF and only keep AF <= 0.05\n maf_df$EAS_AF[is.na(maf_df$EAS_AF)] <- 0\n maf_df <- subset(maf_df,EAS_AF <= 0.05)\n\n # subset Tumor and normal total coverage >=20 & <=1000\n maf_df <- subset(maf_df,t_depth >= 20 & t_depth <= 1000)\n maf_df <- subset(maf_df,n_depth >= 20 & n_depth <= 1000)\n\n # subset Tumor ALT read count >= 3\n maf_df <- subset(maf_df,t_alt_count >= 3)\n\n # subset Normal ALT read count <= 1\n maf_df <- subset(maf_df,n_alt_count <= 1)\n\n # subset Tumor VAF >=0.05\n maf_df$tVAF <- maf_df$t_alt_count\/maf_df$t_depth\n maf_df <- subset(maf_df,tVAF>=0.05)\n\n # subset Normal VAF <=0.05\n maf_df$nVAF <- maf_df$n_alt_count\/maf_df$n_depth\n maf_df <- subset(maf_df,nVAF<=0.05)\n\n maf_df$mut_id <- paste0(maf_df$Tumor_Sample_Barcode,":",\n maf_df$Chromosome,":",\n maf_df$Start_Position,":",\n maf_df$End_Position,":",\n maf_df$Reference_Allele,":",\n maf_df$Tumor_Seq_Allele2)\n\n return(maf_df)\n}\n\nmutect2_filtered <- filter_maf_df(mutect2_maf_df)\nvarscan2_filtered <- filter_maf_df(varscan2_maf_df)\n\nmutect2_varscan2_filtered_intersect <- subset(mutect2_filtered,mut_id %in% varscan2_filtered$mut_id)\ntable(mutect2_varscan2_filtered_intersect$Tumor_Sample_Barcode)\n\nwrite.table(mutect2_filtered,"\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/Mutation\/mutect2_TNpaired\/6patients_mutect2_pass.Somatic.hc.hg38_multianno_filtered.maf",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(varscan2_filtered,"\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/Mutation\/varscan2\/6patients_varscan2_pass.Somatic.hc.hg38_multianno_filtered.maf",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(mutect2_varscan2_filtered_intersect,"\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/Mutation\/6patients_mutect2_varscan2_intersect_pass.Somatic.hc.hg38_multianno_filtered.maf",col.names = T,row.names = F,sep = '\\t',quote = F)\n<\/code><\/pre>\nMerged_maf_df$Tumor_Sample_Barcode <- gsub("-BNPT","",Merged_maf_df$Tumor_Sample_Barcode)\nMerged_maf_df$Tumor_Sample_Barcode <- gsub("-FA","",Merged_maf_df$Tumor_Sample_Barcode)\ntable(Merged_maf_df$Tumor_Sample_Barcode)\n\nwrite.table(Merged_maf_df,"\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/Mutation\/6patients_mutect2_varscan2_intersect_pass_Tumor_Sample_Barcode_merged.Somatic.hc.hg38_multianno_filtered.maf",col.names = T,row.names = F,sep = '\\t',quote = F)<\/code><\/pre>\nmaf2vcf.pl --input-maf 6patients_mutect2_varscan2_intersect_pass_Tumor_Sample_Barcode_merged.Somatic.hc.hg38_multianno_filtered.maf --output-dir \/groups\/phyllodes\/home\/share\/Results\/WES\/maf2vcf --output-vcf 6patients_mutect2_varscan2_intersect_pass_Tumor_Sample_Barcode_merged --ref-fasta \/home\/zhoukaiwen\/database\/GRCh38\/genecode_GRCh38.p14.genome.fa --per-tn-vcfs<\/code><\/pre>\nFETB01-BNPT-E FETB01-B FETB01 FETB01-E\nFETB01-BNPT-M FETB01-B FETB01 FETB01-M\nFETB01-FA-E FETB01-B FETB01 FETB01-E\nFETB01-FA-M FETB01-B FETB01 FETB01-M\nFETB02-BNPT-E FETB02-B FETB02 FETB02-E\nFETB02-BNPT-M FETB02-B FETB02 FETB02-M\nFETB02-FA-E FETB02-B FETB02 FETB02-E\nFETB02-FA-M FETB02-B FETB02 FETB02-M\nFETB03-BNPT-E FETB03-B FETB03 FETB03-E\nFETB03-BNPT-M FETB03-B FETB03 FETB03-M\nFETB03-FA-E FETB03-B FETB03 FETB03-E\nFETB03-FA-M FETB03-B FETB03 FETB03-M\nFETB04-BNPT-E FETB04-B FETB04 FETB04-E\nFETB04-BNPT-M FETB04-B FETB04 FETB04-M\nFETB04-FA-E FETB04-B FETB04 FETB04-E\nFETB04-FA-M FETB04-B FETB04 FETB04-M\nFETB05-BNPT-E FETB05-B FETB05 FETB05-E\nFETB05-BNPT-M FETB05-B FETB05 FETB05-M\nFETB05-FA-E FETB05-B FETB05 FETB05-E\nFETB05-FA-M FETB05-B FETB05 FETB05-M\nFETB06-BNPT-E FETB06-B FETB06 FETB06-E\nFETB06-BNPT-M FETB06-B FETB06 FETB06-M\nFETB06-FA-E FETB06-B FETB06 FETB06-E\nFETB06-FA-M FETB06-B FETB06 FETB06-M<\/code><\/pre>\n# Call \u7a81\u53d8\nperl -ne 'chomp; next if \/^$\/; @a = split \/\\t\/; print "bcftools sort $a[3]\\_vs\\_$a[1].vcf -o $a[3]\\_vs\\_$a[1].sorted.vcf && gatk IndexFeatureFile -I $a[3]\\_vs\\_$a[1].sorted.vcf && gatk Mutect2 -R ~\/database\/GRCh38\/genecode_GRCh38.p14.genome.fa -I ..\/align\/$a[0]_bqsr.bam -I ..\/align\/$a[1]_bqsr.bam -tumor $a[0] -normal $a[1] --alleles $a[3]\\_vs\\_$a[1].sorted.vcf -L $a[3]\\_vs\\_$a[1].sorted.vcf --germline-resource ~\/database\/gatk_bundle\/hg38\/somatic-hg38_af-only-gnomad.hg38.vcf.gz -pon ~\/database\/gatk_bundle\/hg38\/somatic-hg38_1000g_pon.hg38.vcf.gz --f1r2-tar-gz $a[0].shared.sites.f1r2.tar.gz -O $a[0].shared.sites.raw.vcf && echo $a[0] Mutect2 force call ok\\n"' sample_pair.txt > Mutect2_force_call_short.sh\n\n# Annovar\u3001vcf2maf\u7701\u7565\n\n# \u5408\u5e76maf\u6587\u4ef6\ncat *.shared.sites.maf | awk 'NR<=3 || (!\/^#\/ && $1!="Hugo_Symbol")' > 6patients.shared.sites.maf<\/code><\/pre>\n
\n\u63a5\u4e0b\u6765\uff0c\u5236\u4f5cCONIPHER\u8f93\u5165\u6587\u4ef6<\/p>\n# R\nmutation_maf <- read.table("6patients.shared.sites.maf",header = T, sep = '\\t',quote = "")\n\nmutation_maf <- as_tibble(mutation_maf) %>%\n separate(Tumor_Sample_Barcode, into = c("SampleID", "middle", "Tissue"), sep = "-", fill = "right", remove = FALSE) %>%\n mutate(Tissue = coalesce(Tissue, middle)) %>%\n dplyr::select(-middle)\n\nFiltered_mutationTable <- mutation_maf[,c("Tumor_Sample_Barcode","SampleID","Tissue","Chromosome","Start_Position","Reference_Allele","Tumor_Seq_Allele2","t_ref_count","t_alt_count","t_depth")]\n\ncolnames(Filtered_mutationTable) <- c("SAMPLE","CASE_ID","Tissue","CHR","POS","REF","ALT","REF_COUNT","VAR_COUNT","DEPTH")\n\nFiltered_mutationTable$CHR <- gsub("chr","",Filtered_mutationTable$CHR)\n\n# \u6b64\u5904\u9700\u8981\u8fd0\u884csequenza\u4e4b\u540e\u5408\u5e76\u7ed3\u679c\uff0c\u8be6\u89c1sequenza\u6d41\u7a0b\u6587\u4ef6\nmerged_pp <- read.table("\/data02\/zhangmengmeng\/BPT\/LCM-WES\/sequenza\/merged_solutions.txt",header = T,sep = '\\t')\nmerged_seg <- read.table("\/data02\/zhangmengmeng\/BPT\/LCM-WES\/sequenza\/merged_segments.txt",header = T,sep = '\\t')\nmerged_pp_seg <- merge(merged_seg, merged_pp, by = "sample_id", all.x = TRUE)\n\nmerged_pp_seg <- merged_pp_seg[c("sample_id","chromosome","start.pos","end.pos","A","B","cellularity","ploidy")]\nmerged_pp_seg$chromosome <- gsub("chr","",merged_pp_seg$chromosome)\n\nmatched <- Filtered_mutationTable %>%\n inner_join(merged_pp_seg, by = c("SAMPLE" = "sample_id")) %>%\n filter(CHR == chromosome & POS >= start.pos & POS <= end.pos)\n\nmatched <- matched[,-(11:13)]\n\ncolnames(matched)[which(colnames(matched)=="A")] <- "COPY_NUMBER_A"\ncolnames(matched)[which(colnames(matched)=="B")] <- "COPY_NUMBER_B"\ncolnames(matched)[which(colnames(matched)=="cellularity")] <- "ACF"\ncolnames(matched)[which(colnames(matched)=="ploidy")] <- "PLOIDY"\n\n# Then, do an anti_join to get rows from mutationTable that did NOT match\nunmatched <- anti_join(Filtered_mutationTable, matched, by = colnames(Filtered_mutationTable)[1:10])\nunmatched$COPY_NUMBER_A <- 1\nunmatched$COPY_NUMBER_B <- 1\nunmatched$ACF <- "UnKnown"\nunmatched$PLOIDY <- "UnKnown"\n\n# Merge Final input fo CONIPHER #\ntable(colnames(matched) %in% colnames(unmatched))\nConipher_Input <- rbind(matched,unmatched)\nConipher_Input <- Conipher_Input[order(Conipher_Input$SAMPLE),]\nConipher_Input <- Conipher_Input %>%\n group_by(SAMPLE) %>%\n mutate(\n ACF = if_else(ACF == "UnKnown", ACF[ACF != "UnKnown"][1], ACF),\n PLOIDY = if_else(PLOIDY == "UnKnown", PLOIDY[PLOIDY != "UnKnown"][1], PLOIDY)\n ) %>%\n ungroup()\n\nConipher_Input <- Conipher_Input[, c(2, 1, 3:ncol(Conipher_Input))]\nConipher_Input <- subset(Conipher_Input,!(SAMPLE %in% c("FETB01-BPT-M-Unstained","FETB01-BPT-E-Unstained")))\nConipher_Input$SAMPLE <- gsub("-","_",Conipher_Input$SAMPLE)\nConipher_Input$CHR <- factor(Conipher_Input$CHR,levels=c(1:22,"X"))\nConipher_Input <- Conipher_Input[order(Conipher_Input$SAMPLE,Conipher_Input$CHR),]\nConipher_Input <- unique(Conipher_Input)\n\nConipher_Input_FETB01_Epi <- subset(Conipher_Input,CASE_ID == "FETB01" & Tissue =="E")\nConipher_Input_FETB01_Stroma <- subset(Conipher_Input,CASE_ID == "FETB01" & Tissue =="M")\nConipher_Input_FETB02_Epi <- subset(Conipher_Input,CASE_ID == "FETB02" & Tissue =="E")\nConipher_Input_FETB02_Stroma <- subset(Conipher_Input,CASE_ID == "FETB02" & Tissue =="M")\nConipher_Input_FETB03_Epi <- subset(Conipher_Input,CASE_ID == "FETB03" & Tissue =="E")\nConipher_Input_FETB03_Stroma <- subset(Conipher_Input,CASE_ID == "FETB03" & Tissue =="M")\nConipher_Input_FETB04_Epi <- subset(Conipher_Input,CASE_ID == "FETB04" & Tissue =="E")\nConipher_Input_FETB04_Stroma <- subset(Conipher_Input,CASE_ID == "FETB04" & Tissue =="M")\nConipher_Input_FETB05_Epi <- subset(Conipher_Input,CASE_ID == "FETB05" & Tissue =="E")\nConipher_Input_FETB05_Stroma <- subset(Conipher_Input,CASE_ID == "FETB05" & Tissue =="M")\nConipher_Input_FETB06_Epi <- subset(Conipher_Input,CASE_ID == "FETB06" & Tissue =="E")\nConipher_Input_FETB06_Stroma <- subset(Conipher_Input,CASE_ID == "FETB06" & Tissue =="M")\n\nwrite.table(Conipher_Input_FETB01_Epi,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB01_Epi.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB01_Stroma,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB01_Stroma.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB02_Epi,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB02_Epi.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB02_Stroma,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB02_Stroma.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB03_Epi,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB03_Epi.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB03_Stroma,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB03_Stroma.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB04_Epi,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB04_Epi.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB04_Stroma,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB04_Stroma.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB05_Epi,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB05_Epi.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB05_Stroma,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB05_Stroma.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB06_Epi,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB06_Epi.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nwrite.table(Conipher_Input_FETB06_Stroma,"\/data02\/zhangmengmeng\/BPT\/LCM-WES\/conipher\/Conipher_Input_FETB06_Stroma.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\n<\/code><\/pre>\n2.2 \u6b63\u5f0f\u8fd0\u884c<\/h3>\n
# \u6a21\u677f\u5982\u4e0b\nlibrary(CONIPHER)\nlibrary(tidyverse)\n\n# cluster+tree-building\nconipher_run(case_id = "FETB01",\n prefix = "FETB01",\n out_dir = "FETB01",\n input_tsv_loc = "Conipher_Input_FETB01.txt")\n<\/code><\/pre>\n3. \u8fd0\u884c ASCAT<\/h2>\n
# EXAMPLE_ASCAT.R\nlibrary(ASCAT)\noptions(bitmapType='cairo')\nSys.setenv(LD_LIBRARY_PATH = paste("\/data02\/zhangmengmeng\/software\/alleleCount-4.2.1\/lib", Sys.getenv("LD_LIBRARY_PATH"), sep=":"))\nSys.getenv("LD_LIBRARY_PATH")\ndyn.load("\/data02\/zhangmengmeng\/software\/alleleCount-4.2.1\/lib\/libhts.so.3")\nascat.prepareHTS(\n tumourseqfile = "\/data02\/zhangmengmeng\/BPT\/LCM-WES\/align\/TUMOR_bqsr.bam",\n normalseqfile = "\/data02\/zhangmengmeng\/BPT\/LCM-WES\/align\/NORMAL_bqsr.bam",\n tumourname = "TUMOR",\n normalname = "NORMAL",\n allelecounter_exe = "\/data02\/zhangmengmeng\/software\/alleleCount-4.2.1\/bin\/alleleCounter",\n alleles.prefix = "\/data02\/zhangmengmeng\/database\/ASCAT_ref\/G1000_allelesAll_hg38\/G1000_alleles_hg38_chr",\n loci.prefix = "\/data02\/zhangmengmeng\/database\/ASCAT_ref\/G1000_lociAll_hg38\/G1000_loci_hg38_chr",\n BED_file = "\/data02\/zhangmengmeng\/database\/gatk_resource_bundle\/hg38\/AgilentV6_GRCh38_ex_region.sort.bed",\n gender = "XX",\n genomeVersion = "hg38",\n nthreads = 8,\n tumourLogR_file = "TUMOR_LogR.txt",\n tumourBAF_file = "TUMOR_BAF.txt",\n normalLogR_file = "NORMAL_LogR.txt",\n normalBAF_file = "NORMAL_BAF.txt")\n\nascat.bc = ascat.loadData(Tumor_LogR_file = "TUMOR_LogR.txt", Tumor_BAF_file = "TUMOR_BAF.txt", Germline_LogR_file = "NORMAL_LogR.txt", Germline_BAF_file = "NORMAL_BAF.txt", gender = "XX", genomeVersion = "hg38")\nascat.plotRawData(ascat.bc, img.prefix = "TUMOR_Before_correction_")\nascat.bc = ascat.correctLogR(ascat.bc, GCcontentfile = "\/data02\/zhangmengmeng\/database\/ASCAT_ref\/GC_G1000_hg38.txt", replictimingfile = "\/data02\/zhangmengmeng\/database\/ASCAT_ref\/RT_G1000_hg38.txt")\nascat.plotRawData(ascat.bc, img.prefix = "TUMOR_After_correction_")\nascat.bc = ascat.aspcf(ascat.bc)\nascat.plotSegmentedData(ascat.bc)\nascat.output = ascat.runAscat(ascat.bc, gamma=1, write_segments = TRUE)\n\n# Function to count SNPs within segment boundaries\ncount_snps_in_segment <- function(chr, start, end, snp_df) {\n sum(snp_df$Chromosome == chr & snp_df$Position >= start & snp_df$Position <= end)\n}\n\n# Apply function to each row in segments dataframe\nascat.output$segments$nProbes <- mapply(\n count_snps_in_segment,\n ascat.output$segments$chr,\n ascat.output$segments$startpos,\n ascat.output$segments$endpos,\n MoreArgs = list(snp_df = ascat.bc$SNPpo)\n)\nfinal_segment <- ascat.output$segments\ncolnames(final_segment) <- c("SampleID","Chr","Start","End","nMajor","nMinor","nProbes")\nwrite.table(final_segment,"TUMOR_final_segment.txt",col.names = T,row.names = F,sep = '\\t',quote = F)\nQC = ascat.metrics(ascat.bc,ascat.output)\nsave(ascat.bc, ascat.output, QC, file = 'TUMOR_ASCAT_objects.Rdata')<\/code><\/pre>\n# Create_Run_ASCAT.sh\n#!\/bin\/bash\n\n# Path to the example R script\nEXAMPLE_SCRIPT="EXAMPLE_ASCAT.R"\nOUTPUT_SCRIPT="run_ascat.sh"\n\n# Start the script file\necho "#!\/bin\/bash" > "$OUTPUT_SCRIPT"\n\n# Debugging: Check if sample pairs are correctly read\necho "Processing sample pairs..."\n\n# Ensure strict tab-separated parsing\nwhile IFS=$'\\t' read -r TUMOR NORMAL PATIENT_ID; do\n # Trim spaces in TUMOR and NORMAL\n TUMOR=$(echo "$TUMOR" | xargs)\n NORMAL=$(echo "$NORMAL" | xargs)\n PATIENT_ID=$(echo "$PATIENT_ID" | xargs)\n\n # Debug: Show extracted values\n echo "TUMOR: [$TUMOR], NORMAL: [$NORMAL], PATIENT_ID: [$PATIENT_ID]"\n\n # Skip invalid lines\n if [[ -z "$TUMOR" || -z "$NORMAL" ]]; then\n echo "Skipping invalid line: [$TUMOR] [$NORMAL]"\n continue\n fi\n\n # Define output R script name\n R_SCRIPT="ASCAT_${TUMOR}.R"\n\n # Swap TUMOR and NORMAL in the script and create a new file\n awk -v t="$TUMOR" -v n="$NORMAL" '{gsub("TUMOR", t); gsub("NORMAL", n)} 1' "$EXAMPLE_SCRIPT" > "$R_SCRIPT"\n\n # Append the command to run the script (but don't execute)\n echo "\/opt\/R\/4.3.0\/lib64\/R\/bin\/Rscript $R_SCRIPT" >> "$OUTPUT_SCRIPT"\n\ndone < sample_pair.txt\n\n# Make the script executable\nchmod +x "$OUTPUT_SCRIPT"\n\necho "Commands saved in run_ascat.sh. You can run it manually with: bash run_ascat.sh"<\/code><\/pre>\n4. \u8fd0\u884crefphase<\/h2>\n
# R\nlibrary(ASCAT)\nlibrary(refphase)\noptions(bitmapType='cairo')\npatient <- "FETB01-M"\nnormal_sample <- "FETB01-B"\ntumor_samples <- c("FETB01-BNPT-M", "FETB01-FA-M")\n\n# Do some extra work to create a table that will be needed when we run refphase in the next step. This is not needed by ASCAT itself.\nrefphase_sample_data <- data.frame(sample_id = tumor_samples,\n segmentation = paste0(tumor_samples, "_refphase_segs.tsv"),\n snps = paste0(tumor_samples, "_refphase_snps.tsv"),\n purity = NA, ploidy = NA, row.names = tumor_samples)\n\n# ascat_input and ascat_output are required for refphase later on\nascat_input <- list()\nascat_output <- list()\nfor (tumor_sample in tumor_samples) {\n cur_ascat_input <- ascat.loadData(Tumor_LogR_file = paste0("\/data02\/zhangmengmeng\/BPT\/LCM-WES\/ASCAT\/",tumor_sample, "_LogR.txt"),\n Tumor_BAF_file = paste0("\/data02\/zhangmengmeng\/BPT\/LCM-WES\/ASCAT\/",tumor_sample, "_BAF.txt"),\n Germline_LogR_file = paste0("\/data02\/zhangmengmeng\/BPT\/LCM-WES\/ASCAT\/",tumor_sample, "_LogR.txt"),\n Germline_BAF_file = paste0("\/data02\/zhangmengmeng\/BPT\/LCM-WES\/ASCAT\/",tumor_sample, "_BAF.txt"),\n gender = "XX", genomeVersion = "hg38")\n ascat.plotRawData(cur_ascat_input,img.prefix =paste0(tumor_sample,"_Before_correction_"))\n cur_ascat_input = ascat.correctLogR(cur_ascat_input,\n GCcontentfile = "\/data02\/zhangmengmeng\/database\/ASCAT_ref\/GC_G1000_hg38.txt",\n replictimingfile = "\/data02\/zhangmengmeng\/database\/ASCAT_ref\/RT_G1000_hg38.txt")\n ascat.plotRawData(cur_ascat_input,img.prefix =paste0(tumor_sample,"_After_correction_"))\n cur_ascat_input <- ascat.aspcf(cur_ascat_input)\n ascat.plotSegmentedData(cur_ascat_input)\n cur_ascat_output <- ascat.runAscat(cur_ascat_input, gamma = 1.0)\n\n ascat_input[[tumor_sample]] <- cur_ascat_input\n ascat_output[[tumor_sample]] <- cur_ascat_output\n\n # Save the segmentation in the default ASCAT format\n write.table(cur_ascat_output$segments, file = paste0(tumor_sample, "_segments.txt"), sep = "\\t", quote = FALSE, row.names = FALSE)\n\n refphase_sample_data[tumor_sample, "purity"] <- cur_ascat_output$aberrantcellfraction[[1]]\n refphase_sample_data[tumor_sample, "ploidy"] <- cur_ascat_output$ploidy[[1]]\n}\n\nwrite.table(refphase_sample_data, file = paste0(patient,"_refphase-sample-data.tsv"), sep = "\\t", quote = FALSE, row.names = FALSE)\n\n# Load the ascat input and output\nrefphase_input <- refphase_load(\n data_format = "ascat3", samples = tumor_samples,\n ascat_input = ascat_input, ascat_output = ascat_output\n)\n# (Optional) If your data shows reference bias, which presents itself as BAFs\n# in regions with a balanced copy number that are systematically shifted away\n# from 0.5 (usually something like 0.47), this can try to correct for that.\nrefphase_input <- center_baf(refphase_input)\n\n# (Optional) Fit SNP logr data to improve copy number re-estimation in refphase,\n# when using the default ASCAT formula-based method fo re-estimating copy\n# numbers\nrefphase_input <- fit_logr_to_ascat(refphase_input)\n\n# Run refphase on the experimental data\nresults <- refphase(refphase_input)\nsave(results, file = paste0(patient, "-refphase-results.RData"))\n\nwrite_segs(results$phased_segs, file = paste0(patient, "-refphase-segmentation.tsv"))\n\n# (optional) output the SNPs, including phasing information\nwrite_snps(results$phased_snps, file = paste0(patient, "-refphase-phased-snps.tsv.gz"))\n\n# Ploidy might have changed, if we have updated any copy numbers\nwrite.table(results$sample_data, file = paste0(patient, "-refphase-sample-data-updated.tsv"), sep = "\\t", row.names = FALSE)\n\n# Plot the results as PDFs (these files can be huge, over 100 MB in some cases)\npdf(paste0(patient, "-refphase-genome.pdf"), width = 10, height = 4 * nrow(results$sample_data), family = "sans")\nplot(results, what = "genome")\ndev.off()\n\npdf(paste0(patient, "-refphase-chromosomes.pdf"), width = 10, height = 2 + 4 * nrow(results$sample_data), family = "sans")\nplot_all_chromosomes(results$sample_data, results$phased_snps, list(refphase = results$phased_segs, orig = refphase_input$segs), cn_event_calls = results$cn_event_calls)\ndev.off()\n\n# Plot the results as individual PNG files\npng(paste0(patient, "-refphase-genome.png"), width = 600, height = 200 + 300 * nrow(results$sample_data), family = "sans")\nplot(results, what = "genome")\ndev.off()\n\npng(paste0(patient, "-refphase-chromosomes-%02d.png"), width = 700, height = 300 * nrow(results$sample_data), family = "sans")\nplot_all_chromosomes(results$sample_data, results$phased_snps, list(refphase = results$phased_segs, orig = refphase_input$segs), cn_event_calls = results$cn_event_calls)\ndev.off()<\/code><\/pre>\n5. \u8fd0\u884calpaca<\/h2>\n
# R \u53c2\u8003\u6a21\u677f\nload("\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/refphase\/FETB04-M-refphase-results.RData")\ntree_rds <- readRDS("\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/CONIPHER\/FETB04_S\/Trees\/FETB04-M.tree.RDS")\n\nrename_conipher_output <- function(rds){\n rownames(rds$ccf_table_pyclone) <- sub("^([^:]+):", "\\\\1-M:", rownames(rds$ccf_table_pyclone))\n colnames(rds$ccf_table_pyclone) <- sub("_", "-", colnames(rds$ccf_table_pyclone))\n colnames(rds$ccf_table_pyclone) <- sub("_", "-", colnames(rds$ccf_table_pyclone))\n\n rownames(rds$ccf_table_absolute) <- sub("^([^:]+):", "\\\\1-M:", rownames(rds$ccf_table_absolute))\n colnames(rds$ccf_table_absolute) <- sub("_", "-", colnames(rds$ccf_table_absolute))\n colnames(rds$ccf_table_absolute) <- sub("_", "-", colnames(rds$ccf_table_absolute))\n\n rownames(rds$ccf_table_pyclone_clean) <- sub("^([^:]+):", "\\\\1-M:", rownames(rds$ccf_table_pyclone_clean))\n colnames(rds$ccf_table_pyclone_clean) <- sub("_", "-", colnames(rds$ccf_table_pyclone_clean))\n colnames(rds$ccf_table_pyclone_clean) <- sub("_", "-", colnames(rds$ccf_table_pyclone_clean)) \n\n rownames(rds$ccf_table_absolute_clean) <- sub("^([^:]+):", "\\\\1-M:", rownames(rds$ccf_table_absolute_clean))\n colnames(rds$ccf_table_absolute_clean) <- sub("_", "-", colnames(rds$ccf_table_absolute_clean))\n colnames(rds$ccf_table_absolute_clean) <- sub("_", "-", colnames(rds$ccf_table_absolute_clean)) \n\n colnames(rds$clonality_table) <- gsub("_", "-", colnames(rds$clonality_table)) \n\n rds$graph_pyclone$sampleID <- paste0(rds$graph_pyclone$sampleID,"-M")\n rds$parameters$sampleID <- paste0(rds$parameters$sampleID,"-M")\n rds$parameters$prefix <- paste0(rds$parameters$prefix,"-M")\n\n colnames(rds$clonality_out$clonality_table_corrected) <- gsub("_", "-", colnames(rds$clonality_out$clonality_table_corrected)) \n colnames(rds$clonality_out$clonality_table_original) <- gsub("_", "-", colnames(rds$clonality_out$clonality_table_original)) \n colnames(rds$clone_proportion_out$clone_proportion_table) <- gsub("_", "-", colnames(rds$clone_proportion_out$clone_proportion_table)) \n\n rds$subclonal_expansion_score_out$subclonal_exp_score$sample <- gsub("_", "-", rds$subclonal_expansion_score_out$subclonal_exp_score$sample)\n rownames(rds$subclonal_expansion_score_out$subclonal_exp_score) <- rds$subclonal_expansion_score_out$subclonal_exp_score$sample\n\n rds$subclonal_expansion_score_out$subclonal_exp_score_min_sce_trees$`1`$sample <- gsub("_", "-", rds$subclonal_expansion_score_out$subclonal_exp_score_min_sce_trees$`1`$sample)\n rownames(rds$subclonal_expansion_score_out$subclonal_exp_score_min_sce_trees$`1`) <- rds$subclonal_expansion_score_out$subclonal_exp_score_min_sce_trees$`1`$sample\n return(rds)\n}\n\ntree_rds_renamed <- rename_conipher_output(tree_rds)\n\nsaveRDS(tree_rds_renamed,"\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/ALPACA\/FETB04-M.renamed.tree.RDS")<\/code><\/pre>\ntumour_id="FETB04-M"\nrefphase_rData="\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/refphase\/${tumour_id}-refphase-results.RData"\nCONIPHER_tree_object="\/Users\/zhoukaiwen\/Desktop\/Breast_Phyllodes_Tumor\/Bioinfo\/WES\/ALPACA\/${tumour_id}.renamed.tree.RDS"\nprep_input_dir=".\/${tumour_id}\/input"\nrun_input_dir="${tumour_id}\/input"\nrun_output_dir="${tumour_id}\/output"\n\nalpaca input-conversion \\\n --tumour_id $tumour_id \\\n --refphase_rData $refphase_rData \\\n --CONIPHER_tree_object $CONIPHER_tree_object \\\n --CONIPHER_tree_index 1 \\\n --output_dir $prep_input_dir\n\n# \u91cd\u65b0\u8c03\u6574\u6837\u672c\u547d\u540d\uff0cR\u8f93\u51fa\u65f6\u81ea\u52a8\u4f1a\u628a\u5217\u540d\u4e2d\u7684-\u4fee\u6539\u4e3a_\uff0c\u8fd9\u91cc\u6539\u56de\u6765\nsed -i '' 's\/_\/-\/g' "$(pwd)\/${run_input_dir}\/cp_table.csv"\n\n# \u4f1a\u6709\u4e0d\u91cd\u590d\u7684CN segments\uff0c\u8fd9\u91cc\u53ea\u6311\u9009\u91cd\u590d\u7684CNV\u533a\u6bb5\nmv "$(pwd)\/${run_input_dir}\/ALPACA_input_table.csv" "$(pwd)\/${run_input_dir}\/ALPACA_input_table.csv.backup"\nawk -F, 'NR==1{print; next} {count[$3]++; lines[NR]=$0; seg[NR]=$3} END{for(i=2;i<=NR;i++) if(count[seg[i]]>1) print lines[i]}' "$(pwd)\/${run_input_dir}\/ALPACA_input_table.csv.backup" > "$(pwd)\/${run_input_dir}\/ALPACA_input_table.csv"\n\n# \u8f6f\u4ef6bug\uff1aalpaca run\u8bbe\u7f6e\u4e86run_input_dir\u7684\u60c5\u51b5\u4e0b\uff0c\u4f1a\u9ed8\u8ba4\u8f93\u5165\u8def\u5f84\u4e3a${tumour_id}\/${tumour_id}\uff0c\u6545\u8fd0\u884c\u524d\u521b\u5efa\u8f6f\u94fe\u63a5\nln -s "$(pwd)\/${tumour_id}\/input" "${tumour_id}\/${tumour_id}"\n\nalpaca run \\\n --input_tumour_directory "${run_input_dir}" \\\n --output_directory "${run_output_dir}"<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"