Call cnvs using cn.mops #cluster

cn.mops.R

#!/usr/bin/Rscript

#SBATCH --job-name=mops

#SBATCH --output=../log/%j.txt
#SBATCH --error=../log/%j.out

#SBATCH --partition=compute

#SBATCH --nodelist=node3
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=16
#SBATCH --nodes=1
#SBATCH --mem=16384

#SBATCH [email protected]
#SBATCH --workdir=/lscr2/andersenlab/dec211/RUN/v2_snpset/bam/

# This script will call copy number variants using cn.mops

library(cn.mops)
library(dplyr)

BAMFiles <- list.files(pattern=".bam$")
bamDataRanges <- getReadCountsFromBAM(BAMFiles, mode="paired", refSeqName=c("I","II","III","IV","V","X","MtDNA"), parallel=12)
results <- calcIntegerCopyNumbers(cn.mops(bamDataRanges, parallel=16, returnPosterior=TRUE))
results_haplo <- calcIntegerCopyNumbers(haplocn.mops(bamDataRanges, parallel=16))

# Plot results for each sample
for (i in 1:length(BAMFiles)) {
	png(paste0("../indel/cn.mops/",BAMFiles[i], ".png"), width=1800, height=1200)
		segplot(results, sampleIdx=i)
	dev.off()
}

for (i in 1:length(BAMFiles)) {
	png(paste0("../indel/cn.mops/haplo_",BAMFiles[i], "_haplo.png"), width=1800, height=1200)
		segplot(results_haplo, sampleIdx=i)
	dev.off()
}

# Save Rdata
save(results, "../indel/cn.mops/results.Rdata")
save(results_haplo, "../indel/cn.mops/results_haplo.Rdata")
save(bamDataRanges, "../indel/cn.mops/bamDataRanges.Rdata")

fix_name <- function(x) { gsub("CN","",x) }

# Diploid Results
results_df <- tbl_df(as.data.frame(cnvr(results))) %>%
				   rename(chrom=seqnames) %>%
				   mutate(name="", score=0) %>%
				   mutate_each(funs(fix_name), matches("bam")) %>%
				   select(chrom, start, end,  name, score, strand, width, everything()) %>%
				   arrange(chrom) 
# Fix names
names(results_df) <- gsub("\\.bam","",names(results_df))
write.table(results_df, file="../indel/cn_mops_df.bed", quote=F, sep="\t", row.names=F, col.names=T)

# Haploid Results

results_haplo_df <- tbl_df(as.data.frame(cnvr(results_haplo))) %>%
				   rename(chrom=seqnames) %>%
				   mutate(name="", score=0) %>%
				   mutate_each(funs(fix_name), matches("bam")) %>%
				   select(chrom, start, end,  name, score, strand, width, everything()) %>%
				   arrange(chrom) 
# Fix names
names(results_haplo_df) <- gsub("\\.bam","",names(results_haplo_df))
write.table(results_haplo_df, file="../indel/cn_mops_haplo_df.bed", quote=F, sep="\t", row.names=F, col.names=T)