mappabilityFile = "wgEncodeCrgMapabilityAlign100mer.bigWig.Window.500.start.1.end.48129895.txt",
August 11, 2016 ยท View on GitHub
This is a faster way to install/check results of SRBreak
Users need to install samtools if required.
git clone https://github.com/hoangtn/CNVrd2.git
git clone https://github.com/hoangtn/SRBreak.git
cd SRBreak/
cat SRBreak/R/SRBreak.R > ../allSourceFile.R
cd ../CNVrd2/
cat CNVrd2/R/*R >> ../allSourceFile.R
#sed -i .bak "s/CNVRd2:://g" ../allSourceFile.R
#rm ../allSourceFile.R.bak
cd ..
Then, go inside R and install three packages: Rsamtools, VariantAnnotation, BSgenome.Hsapiens.UCSC.hg19, mclust
##Install packages
source("https://bioconductor.org/biocLite.R")
biocLite("Rsamtools")
biocLite("VariantAnnotation")
biocLite("BSgenome.Hsapiens.UCSC.hg19")
biocLite("DNAcopy")
install.packages("mclust")
Test for whole chromosome chr21
Users should download whole-chromosome files from:
https://www.dropbox.com/sh/wejg4r37kdkjokc/AADglfiSzny1bJ4w1ROAVERoa?dl=0
After that, unzip the files into one directory (e.g., Chr21BamFile)
unzip Chr21BamFile.zip -d Chr21BamFile
In R
source("allSourceFile.R")
library("VariantAnnotation")
library("BSgenome.Hsapiens.UCSC.hg19")
library("mclust")
library("Rsamtools")
#link file: https://www.dropbox.com/sh/wejg4r37kdkjokc/AADglfiSzny1bJ4w1ROAVERoa?dl=0
#
ePos <- 48129895 ##End position
windows = 500
dirBamFile = "./Chr21BamFile/"
segmentalDuplicationFile = paste0(dirBamFile, "GRCh37GenomicSuperDup.tab.onlyChr21")
##Segmental duplication information can be obtained from here: http://humanparalogy.gs.washington.edu/build37/data/GRCh37GenomicSuperDup.tab
##and should be formatted as the file GRCh37GenomicSuperDup.tab.onlyChr21 for each chromosome
st = 9500001
en = ePos
st = 9500001 ##Start position; if we don't know, we can set st = 1 and SRBreak can adjust this value inside to remove unknown regions
chr = "chr21" ##Note: these bam files are only from chr21
system.time(outputSRBreak <- SRBreak(readDepthWindow = windows,##read-depth window size
chr = chr, ##Chromosome name
st = st, ##Start position
en = en, ##End position
dirBamFile = dirBamFile, ##Bam files' directory
detectAllRegion = TRUE, ##Set this = TRUE in order to obtain all CNV regions
rdQualityMapping = 0, ##Mapping quality (used in read-depth method)
testType = "Count", ##Test type
correctGC = TRUE, # FALSE, ##Correct GC content
upperCNThreshold = 0.25, ##Larger than this threshold is duplication
lowerCNThreshold = -0.25, ##Smaller than this threshold is deletion
countThreshold = 2, ##Number of duplications/deletions: should be >= 2
minLengthSV = 1000, ##Minimum length of a duplication/deletion event
# mappabilityFile = "wgEncodeCrgMapabilityAlign100mer.bigWig.Window.500.start.1.end.48129895.txt",
##Mappability information is not a compulsory option, it should be turned off
usingPairedEnds = FALSE, ##Not use paired-end information,
thresholdOfIntersectionBetweenRDandPEM = 0.9 #Not useful if only single-end reads used
, segmentalDuplicationFile = segmentalDuplicationFile,
adjustStartPosition = TRUE ##Adjust positions
))
rawoutputSV <- outputSRBreak$svResult
dim(rawoutputSV)
finalResult <- rawoutputSV
dim(finalResult)
This file ./test_chr21_SRBreak.ipynb shows results for the test data. The test results include: whole chr21, a specific region, and a specific region for two samples.