Merge branch 'dev' into peak_calling

bin/bed_to_fasta.R

@@ -0,0 +1,28 @@
+#!/usr/bin/env Rscript
+
+# Splitting BED-files depending on their cluster.
+# The Sequences of each cluster are writen as an FASTA-file.
+# @parameter bedInput <string> BED-file with sequences and cluster-id as column"TEs
+# @parameter prefix <string> prefix for filenames
+# @parameter min_seq <INT> min. number of sequences per cluster
+
+args = commandArgs(trailingOnly = TRUE)
+
+bedInput <- args[1]
+prefix <- args[2] 
+min_seq <- args[3] 
+
+bed <- data.table::fread(bedInput, header = FALSE, sep = "\t")
+
+clusters <- split(bed, bed$V8, sorted = TRUE, flatten = FALSE) # <---- Spalte mit Cluster
+discard <- lapply(1:length(clusters), function(i){
+  clust <- as.data.frame(clusters[i])
+  print(nrow(clust))
+  if (nrow(clust) >= as.numeric(min_seq) ) {
+    sequences <- as.list(clust[[7]])   # <---- Splate mit Sequenz
+    outfile <- paste0(prefix,"_cluster_",i,".FASTA")
+    seqinr::write.fasta(sequences = sequences, names = clust[[4]], file.out = outfile, as.string = TRUE) # <---- Spalte mit Name
+  } else {
+    print(paste0("Cluster: ",i," is to small"))
+  }
+})

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 \<path/to/file.bed\> --motifs \<path/to/known/motifs.bed\> --fasta \<path/to/genome.fasta\> \[optional_parameter \<value\>\] ..."
+	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 \<path/data.bed\> --motifs \<path/motifs.bed\> --fasta \<path/file.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

bin/get_best_motif.py

@@ -0,0 +1,23 @@
+# parses arguments using argparse
+# @return args list of all parameters
+def parse_arguments():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("meme", help="Path to meme file")
+    parser.add_argument("output", help="")
+    args = parser.parse_args()
+    return args
+
+
+def main():
+    args = parse_arguments()
+    out = open(args.output, "w+")
+    with open(args.meme) as f:
+        for line in f:
+            if 'MOTIF  2' in line:
+                break
+            out.write(line)
+
+
+if __name__ == "__main__":
+    import argparse
+    main()

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')

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')
+