Implemented filtering of CNV-variable regions.

MoCaSeq.sh

@@ -871,7 +871,7 @@ if [ $sequencing_type = 'WES' ]; then
 	echo '---- Run CopywriteR ----' | tee -a $name/results/QC/$name.report.txt
 	echo -e "$(date) \t timestamp: $(date +%s)" | tee -a $name/results/QC/$name.report.txt
 
-	Rscript $repository_dir/CNV_RunCopywriter.R $name $species $threads $runmode $genome_dir $types
+	Rscript $repository_dir/CNV_RunCopywriter.R $name $species $threads $runmode $genome_dir $centromere_file $varregions_file $types
 
 	echo '---- Export raw data and re-normalize using Mode ----' | tee -a $name/results/QC/$name.report.txt
 	echo -e "$(date) \t timestamp: $(date +%s)" | tee -a $name/results/QC/$name.report.txt

repository/CNV_RunCopywriter.R

@@ -15,9 +15,14 @@ species = args[2]
 threads = args[3]
 runmode = args[4]
 genome_dir = args[5]
-types =  args[6]
+centromere_file <- args[6]
+varregions_file <- args[7]
+types =  args[8]
 
-library("CopywriteR")
+library(CopywriteR)
+library(GenomeInfoDb)
+library(naturalsort)
+library(GenomicRanges)
 
 tumor_bam = paste(name,"/results/bam/",name,".Tumor.bam",sep="")
 normal_bam = paste(name,"/results/bam/",name,".Normal.bam",sep="")
@@ -47,4 +52,38 @@ CopywriteR(sample.control = sample.control,
              reference.folder = file.path(reference_files),
              bp.param = bp.param)
 
+log2.reads=read.table(paste0(name,"/results/Copywriter/CNAprofiles/log2_read_counts.igv"), header=T, sep="\t",check.names=FALSE)
+
+file.copy(paste0(name,"/results/Copywriter/CNAprofiles/log2_read_counts.igv"),paste0(name,"/results/Copywriter/CNAprofiles/log2_read_counts_backup.igv"), overwrite=T)
+
+log2.reads.GR=GRanges(log2.reads$Chromosome, IRanges(log2.reads$Start, log2.reads$End),Feature=as.character(log2.reads$Feature), Normal=log2.reads[,5],Tumor=log2.reads[,6])
+
+# remove regions with increased variability for mice and centromere regions for humams
+if (species == "Human")
+{
+	filter=read.delim(centromere_file)
+	flankLength=5000000
+}
+if (species == "Mouse")
+{
+	filter=read.delim(varregions_file)
+	flankLength=0
+}
+
+colnames(filter)[1:3] <- c("Chromosome","Start","End")
+filter$Start <- filter$Start - flankLength
+filter$End <- filter$End + flankLength
+filter=GRanges(filter$Chromosome, IRanges(filter$Start, filter$End))
+
+hits <- findOverlaps(query = log2.reads.GR, subject = filter)
+ind <- queryHits(hits)
+message("Removed ", length(ind), " bins near centromeres.")
+log2.reads.GR=(log2.reads.GR[-ind, ])
+
+log2.reads.fixed=as.data.frame(log2.reads.GR)
+log2.reads.fixed=log2.reads.fixed[,c("seqnames", "start", "end", "Feature", "Normal","Tumor")]
+colnames(log2.reads.fixed)=c("Chromosome", "Start", "End", "Feature",paste0("log2.",name,".Normal.bam"),paste0("log2.",name,".Tumor.bam"))
+
+write.table(log2.reads.fixed,paste0(name,"/results/Copywriter/CNAprofiles/log2_read_counts.igv"), sep="\t", quote=F, row.names=F,col.names=T)
+
 plotCNA(destination.folder = file.path(paste(name,"/results/Copywriter/",sep="")))