From 87bbda0c4cd6ccd1d5507f54e480aa975f13813d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ren=C3=A9=20Wiegandt?= Date: Thu, 8 Nov 2018 07:26:59 -0500 Subject: [PATCH] added scripts for part 1.2 --- bin/compareBed.sh | 183 ++++++++++++++++++++++++++++++++++++++++++++++ bin/maxScore.R | 9 +++ bin/merge.R | 23 ++++++ 3 files changed, 215 insertions(+) create mode 100644 bin/compareBed.sh create mode 100644 bin/maxScore.R create mode 100644 bin/merge.R diff --git a/bin/compareBed.sh b/bin/compareBed.sh new file mode 100644 index 0000000..bcf52fb --- /dev/null +++ b/bin/compareBed.sh @@ -0,0 +1,183 @@ +#!/bin/bash +#data +#motifs +#workdir +#fasta +#min +#max +#output +wrong=() +da=false +mo=false +wo=false +fa=false +mi=false +ma=false +ou=false +he=false + +if [ $# -eq 0 ] +then + he=true +fi + +while [[ $# -gt 0 ]] +do +key="$1" + +case $key in + -d|--data) + data="$2" + da=true + shift + shift + ;; + -m|--motifs) + motifs="$2" + mo=true + shift + shift + ;; + -w|--workdir) + workdir="$2" + wo=true + shift + shift + ;; + -f|--fasta) + fasta="$2" + fa=true + shift + shift + ;; + -min|--min) + min=$2 + mi=true + shift + shift + ;; + -max|--max) + max=$2 + ma=true + shift + shift + ;; + -o|--output) + output="$2" + ou=true + shift + shift + ;; + -h|--help) + help=true + he=true + shift + ;; + *) + wrong+=("$1") + shift + ;; +esac +done + +count=${#wrong[@]} +if [ $count -gt 0 ] +then + for i in ${wrong[@]} + do + echo wrong parameter $i + echo call script without parameters for help or call --help + echo exit + done +exit 1 +fi + +if [ $he == true ] +then + echo "This script utilies bedtools to select new footprints from data." + echo "Therefore the data is compared with known footprint motifs." + echo "The output is a new .bed file with added sequence information." + echo "Required arguments are data and motifs, both in bed-format, and the fasta genome sequences." + echo "If a parameter is chosen, a value must be provided or an error will occur." + echo "--------------------" + echo "usage: compareBed.sh --data \ --motifs \ --fasta \ \[optional_parameter \\] ..." + echo "--------------------" + echo " required parameter:" + echo " -d --data the path to the .bed file of the footprints" + echo " -m --motifs the path to the .bed file of the scanned motifs" + echo " -f --fasta the path to the .fasta file of genome" + echo " " + echo " optional parameter:" + echo " -w --workdir the path to directory where temporary files will be created" + echo " default is a new directory that is created in the current directory" + echo " -min --min minimum size of footprints\; default is 10" + echo " -max --max maximum size of footprints\; default is 60" + echo " -o --output output path/file \; default dir is workdir and filename is newFootprints.bed and newFootprints.bed.fasta" + echo " -h --help shows this help message" +exit 0 +fi + +echo selected parameters +echo ------------------- +echo data: $da \(required\) +echo motifs: $mo \(required\) +echo workdir: $wo +echo fasta: $fa \(required\) +echo min: $mi +echo max: $ma +echo output: $ou +echo help: $he + +if [ $da == false ] || [ $mo == false ] || [ $fa == false ] +then + echo required parameters not given. + echo required are: --data \ --motifs \ --fasta \ + exit 1 +fi + +if [ $wo == false ] +then + if [ ! -d workdir1337 ] + then + mkdir workdir1337 + fi + wo=true + workdir="workdir1337" +fi + +if [ $ou == false ] +then + output="newFootprints.bed" + ou=true +fi + +if [ $mi == false ] +then + min=10 + mi=true +fi + +if [ $ma == false ] +then + max=60 + ma=true +fi +#1. first filter. no overlap vs. overlap +bedtools intersect -v -a $data -b $motifs > "$workdir"/pass1Tr.bed +bedtools intersect -wa -a $data -b $motifs > "$workdir"/pass1Fa.bed + +#2. compute absolut maxscore position +Rscript --vanilla maxScore.R "$workdir"/pass1Fa.bed + +#3. subtract overlapping regions for additional motifs +bedtools subtract -a "$workdir"/pass1Fa.bed -b $motifs > "$workdir"/pass2Tr.bed + +#4. remove short/long motivs, make unique ids (relevant for some splitted tfbs from subtract) and handle maxScorePosition +Rscript --vanilla merge.R $min $max "$workdir" + +#5. add fasta sequences to bed and create fasta file +bedtools getfasta -fi $fasta -bed "$workdir"/merged.bed -bedOut > "$workdir"/"$output" +bedtools getfasta -name -fi $fasta -bed "$workdir"/"$output" -fo "$workdir"/"$output".fasta + +#6 clean up +rm "$workdir"/pass1Fa.bed "$workdir"/pass1Tr.bed "$workdir"/pass2Tr.bed "$workdir"/merged.bed \ No newline at end of file diff --git a/bin/maxScore.R b/bin/maxScore.R new file mode 100644 index 0000000..4ea545b --- /dev/null +++ b/bin/maxScore.R @@ -0,0 +1,9 @@ +#!/home/jhamp/.conda/envs/tfbs/bin/Rscript +args = commandArgs(TRUE) +file = args[1] + +tab = read.table(file, header=FALSE) +colnames(tab) = c("chromosome", "start", "stop", "id", "score", "maxpos", "length") +tab$maxpos = tab$start + tab$maxpos + +write.table(tab, file, row.names=FALSE, col.names=FALSE, quote=FALSE, sep='\t') diff --git a/bin/merge.R b/bin/merge.R new file mode 100644 index 0000000..bbe3c9a --- /dev/null +++ b/bin/merge.R @@ -0,0 +1,23 @@ +#!/home/jhamp/.conda/envs/tfbs/bin/Rscript +args=commandArgs(TRUE) +min=as.numeric(args[1]) +max=as.numeric(args[2]) +folder=args[3] +splitted = read.table(paste(folder, "/pass2Tr.bed", sep=''), header=FALSE) +colnames(splitted) = c("chromosome", "start", "stop", "id", "score", "maxpos", "length") +p1 = read.table(paste(folder, "/pass1Tr.bed", sep=''), header=FALSE) +colnames(p1) = c("chromosome", "start", "stop", "id", "score", "maxpos", "length") +p1$maxpos = p1$start + p1$maxpos + +splitted=rbind(splitted, p1) + +splitted=splitted[which(splitted$stop - splitted$start >= min),] +splitted=splitted[which(splitted$stop - splitted$start <= max),] +splitted$id=make.unique(as.character(splitted$id)) +splitted$length=splitted$stop - splitted$start + +splitted=cbind(splitted, containsMaxpos=0) +splitted$containsMaxpos[intersect(which(splitted$start <= splitted$maxpos), which(splitted$stop > splitted$maxpos))] = 1 +splitted$maxpos = splitted$maxpos - splitted$start +write.table(splitted, paste(folder, "/merged.bed", sep=''), row.names=FALSE, col.names=FALSE, quote=FALSE, sep='\t') +