CopywriteR code complete

CopywriteR 은 off-target sequenceing 데이터를 이용하여 CNV(Copy number variation) 를 검출하는 알고리즘 이다.

workflow는 크게 3단계로 나누어 진다.

 

preCopywriteR - 파라미터 설정 ((blacklist, bin region, GC-content, and mappability .bed files, reference genomes)

     output.folder - output folder 경로 설정

     bin.size - bin length (in bp) 설정

     ref.genome - 참조서열 설정 (hg18, hg19, hg38, mm9 or mm10.)

     prefix - prefix 설정 (ex - chromosome notation ,prefix = "chr".)

 

CopywriteR - off-target reads 들을 이용해서 DNA CNV 가 생성

     sample.control - sample 과 control BAM 파일을 하나의 데이터프레임으로 생성

     destination.folder - output 파일 경로

     reference.folder - 잠조서열 경로

     bp.param -  BiocParallelParam ( parallel computing 설정)

     capture.regions.file - capture regions 이 .bed format 으로 나오기 위해서 경로 설정.  read 가 Overlapping 되는 부분                                             reduced.

     keep.intermediary.files - 분석이 끝나고 peak regions 이 보관 되도록 하는 옵션 (Defaults to FALSE)

 

plotCNA - whole-genome and per-chromosome plots 생성

     destination.folder - output 파일 경로

     smoothed - copy number smoothens 단계 (Defaults to TRUE.)

     sample.plot - 모든 sample 을 reference 유무에 따라 모두 plotting.

 

위 plot 이 아래 plot 에 비해 bin 의 분포가 퍼져 있는것을 확인 할 수 있다. (MAD = 0.284) 빨간색 줄이 그 bin 들의 평균치를 나타내느 segment 인데 이것 또한 다양한 분포를 가진다. (해당 샘플의 성별은 남자임으로 chrX 와 chrY 에 해당하는 23 과 24 가 1 copy씩 있는걸로 보임)

 

#####################################
# Rscript onlyCopywriteR_test_220921.R
# check path of working directiroy and input bam files
####################################

getwd()
setwd("/data/skkwon1048/Test/test_analysis_copywriteR_220920/")

#BiocManager::install("CopywriteR")

library(CopywriteR)
library(CopyhelpeR)
getPathHelperFiles("hg19")

data.folder <- tools::file_path_as_absolute(file.path(getwd()))

#generate helper files (blacklist, bin region, GC-content and mappability.bed files)
preCopywriteR(output.folder = file.path(data.folder),
  bin.size = 50000,
  ref.genome = "hg19",
  prefix = "chr")

list.dirs(path = file.path(data.folder), full.names = FALSE)[2]

list.files(path = file.path(data.folder, "hg19_50kb_chr"), full.names = FALSE)

load(file = file.path(data.folder, "hg19_50kb_chr", "blacklist.rda"))

blacklist.grange

load(file = file.path(data.folder, "hg19_50kb_chr", "GC_mappability.rda"))
GC.mappa.grange[1001:1011]

bp.param <- SnowParam(workers = 1, type = "SOCK")
bp.param

#CopywriteR()
#BiocManager::install("SCLCBam")

#library(SCLCbam)
#path <- SCLCBam::getPathBamFolder()

path <- "/data/skkwon1048/Test/test_analysis_copywriteR_220920"
samples <- list.files(path = path, pattern = "bam$", full.names = TRUE)

#setting control
controls <- samples
#controls <- samples[rep(2, length(samples))]

sample.control <- data.frame(samples, controls)

CopywriteR(sample.control = sample.control,
  destination.folder = file.path(data.folder),
  reference.folder = file.path(data.folder, "hg19_50kb_chr"),
  bp.param = bp.param)
  
cat(list.files(path = file.path(data.folder, "CNAprofiles")), sep = "\n")

read.table(file = file.path(data.folder, "CNAprofiles", "read_counts.txt"), header = TRUE)[1001:1006, ]

read.table(file = file.path(data.folder, "CNAprofiles",
  "log2_read_counts.igv"), header = TRUE)[817:822, ]
  
cat(list.files(path = file.path(data.folder, "CNAprofiles", "qc")), sep = "\n")

#plotCNA()
#analyses CopywriteR output using the segmentation algorithm CBS and creates both whole-genome and per-chromosome plots.
plotCNA(destination.folder = file.path(data.folder))

cat(list.files(path = file.path(data.folder, "CNAprofiles")), sep = "\n")

cat(list.files(path = file.path(data.folder, "CNAprofiles", "plots")), sep = "\n")

 

input file .bam 과 output file 경로 설정하면 코드 구현이 가능합니다.

 

아래 링크 참고 바랍니다.

 

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0617

CopywriteR: DNA copy number detection from off-target sequence data - Genome Biology

 

https://github.com/PeeperLab/CopywriteR

GitHub - PeeperLab/CopywriteR: DNA copy number detection from off-target sequence data