Commit 451e6cfc authored by Holger Brandl's avatar Holger Brandl
Browse files

updated utilities for sorted overlap analyses

parent 587dfe0b
......@@ -102,7 +102,7 @@ export -f cs_bowtie_qc
cs_lib_size(){
rm algn_counts.txt
rm -f algn_counts.txt
for bamFile in $(ls $*); do
(
echo "processing $bamFile"
......@@ -132,11 +132,11 @@ export -f cs_lib_size
cs_sync_order_chrominfo2bam(){
if [ $# -ne 2 ]; then
echo "Usage: cs_sync_order_bam2bed <bam_file> <chrom_info>" >&2 ; return;
echo "Usage: cs_sync_order_chrominfo2bam <bam_file> <chrom_info>" >&2 ; return;
fi
echo '
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.26/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.38/R/core_commons.R")
bamFile=commandArgs(T)[1]
chromInfoFile=commandArgs(T)[2]
......@@ -147,9 +147,9 @@ chromInfoFile=commandArgs(T)[2]
#bamFile="/projects/bioinfo/holger/projects/krause_chipseq/reanalysis/../alignments_mmfilt/H3HA_Sphere_mmf.bam";
#chromInfoFile="danRer7.chromInfo"
chrOrderBam <- readLines(pipe(paste("samtools view -H ", bamFile, " | cut -f2 | grep -F SN"))) %>% str_replace(fixed("SN:"), "")
chrOrderBam = readLines(pipe(paste("samtools view -H ", bamFile, " | cut -f2 | grep -F SN"))) %>% str_replace(fixed("SN:"), "")
chromInfo <- read_tsv(chromInfoFile, col_names=F) %>% set_names("chromosome", "chr_length") %>% mutate(chrom=str_replace(chromosome, "chr", ""))
#chromInfo = read_tsv(chromInfoFile, col_names=F) %>% subset(select=1:2) %>% set_names("chromosome", "chr_length") #%>% mutate(chrom=str_replace(chromosome, "chr", ""))
#data.frame(chrom=chrOrderBam) %>%
# left_join(chromInfo) %>%
......@@ -159,36 +159,40 @@ chromInfo <- read_tsv(chromInfoFile, col_names=F) %>% set_names("chromosome", "c
read_tsv(chromInfoFile, col_names=F) %>%
set_names("chromosome", "chr_length") %>%
mutate(chromosome=str_replace(chromosome, "M", "MT") %>% str_replace("chr", "")) %>%
mutate(chromosome=str_replace(chromosome, "chrM", "M") %>% str_replace("chr", "")) %>%
subset(select=1:2) %>%
# mutate(chromosome=str_replace(chromosome, "M", "MT") %>% str_replace("chr", "")) %>%
# mutate(chromosome=str_replace(chromosome, "chrM", "M") %>% str_replace("chr", "")) %>%
filter(chromosome %in% chrOrderBam) %>%
mutate(chromosome=factor(as.character(chromosome), levels=chrOrderBam)) %>% arrange(chromosome) %>%
filter(!is.na(chromosome)) %>%
write.delim(stdout(), header=F, quote=F)
' | Rscript --vanilla - $1 $2 2> /dev/null
}
export -f cs_sync_order_chrominfo2bam
cs_sync_order_bam2bed(){
cs_sync_chrominfo_order_to_bed(){
if [ $# -ne 2 ]; then
echo "Usage: cs_sync_order_bam2bed <bam_file> <bed_file>" >&2 ; return;
echo "Usage: cs_sync_chrominfo_order_to_bed <sorted_chrom_info> <bed_file>" >&2 ; return;
fi
echo '
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.26/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.38/R/core_commons.R")
bamFile=commandArgs(T)[1]
chromInfo=commandArgs(T)[1]
bedFile=commandArgs(T)[2]
options(warn = -1)
#bamFile="/projects/bioinfo/holger/projects/khan_chipseq_h1/alignments_mmfilt/H1M_Dome_mmf.bam"
#bedFile="zv9_ens79_genebodies.bed"
chrOrderBam = read_tsv(chromInfo, col_names=F) %>% set_names("chromosome", "length") %$% chromosome
bedData = read.delim(bedFile, h=F)
chrOrderBam <- readLines(pipe(paste("samtools view -H ", bamFile, " | cut -f2 | grep -F SN"))) %>% str_replace(fixed("SN:"), "")
#read.delim(bedFile, h=F) %>% slice(match(chrOrderBam, as.character(V1))) %>%
read.delim(bedFile, h=F) %>% mutate(V1=factor(as.character(V1), levels=chrOrderBam)) %>% arrange(V1, V2) %>% write.delim(stdout(), header=F, quote=F)
bedData %<>% filter(V1 %in% chrOrderBam) %>% mutate(V1=factor(as.character(V1), levels=chrOrderBam)) %>% arrange(V1, V2)
bedData %>% write.delim(stdout(), header=F, quote=F)
' | Rscript --vanilla - $1 $2 2> /dev/null
}
export -f cs_sync_order_bam2bed
export -f cs_sync_chrominfo_order_to_bed
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment