diff --git a/bin/bed_to_fasta.R b/bin/bed_to_fasta.R index 6767ab5..09e3f4f 100644 --- a/bin/bed_to_fasta.R +++ b/bin/bed_to_fasta.R @@ -4,18 +4,25 @@ # The Sequences of each cluster are writen as an FASTA-file. # @parameter bedInput BED-file with sequences and cluster-id as column"TEs # @parameter prefix prefix for filenames +# @parameter min_seq min. number of sequences per cluster args = commandArgs(trailingOnly = TRUE) bedInput <- args[1] -prefix <- args[2] # "Fasta" +prefix <- args[2] +min_seq <- args[3] bed <- data.table::fread(bedInput, header = FALSE, sep = "\t") -clusters <- split(bed, bed$V3, sorted = TRUE, flatten = FALSE) # <---- Spalte mit Cluster +clusters <- split(bed, bed$V8, sorted = TRUE, flatten = FALSE) # <---- Spalte mit Cluster discard <- lapply(1:length(clusters), function(i){ - sequences <- as.list(as.data.frame(clusters[i])[[2]]) # <---- Splate mit Sequenz - outfile <- paste0(prefix,"_cluster_",i,".FASTA") - seqinr::write.fasta(sequences = sequences, names = as.data.frame(clusters[i])[[1]] - , file.out = outfile, as.string = TRUE) # <---- Spalte mit Name + 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")) + } }) 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') + diff --git a/pipeline.nf b/pipeline.nf index 49adb09..e4d502f 100644 --- a/pipeline.nf +++ b/pipeline.nf @@ -1,8 +1,11 @@ //!/usr/bin/env nextflow Channel.fromPath(params.input).map {it -> [it.simpleName, it]}.set {bigwig_input} -Channel.fromPath(params.genome_fasta).into {fa_overlap; fa_scan} +Channel.fromPath(params.genome_fasta).into {fa_overlap; fa_scan; fa_overlap_2} Channel.fromPath(params.jaspar_db).into {db_for_motivscan; db_for_tomtom} +Channel.fromPath(params.config).set {config} + +params.min_seq = 10 process footprint_extraction { @@ -12,14 +15,15 @@ process footprint_extraction { set name, file (bigWig) from bigwig_input output: - set name, file ('*.bed') into bed_for_overlapp_with_TFBS + set name, file ('*.bed') into bed_for_overlap_with_TFBS script: """ + python ${path_bin}/call_peaks.py --bigwig ${bigWig} --bed ${bed} --output_file ${name}_called_peaks.bed --window_length ${params.window_length} --step ${params.step} --threshold ${params.threshold_footprint_extraction} """ } - +//Abfrage ob ausgeführt werden muss. process extract_known_TFBS { input: @@ -35,16 +39,20 @@ process extract_known_TFBS { } +bed_for_overlap_with_TFBS.combine(known_TFBS_for_overlap).combine(fa_overlap_2).set {for_overlap} + + process overlap_with_known_TFBS { input: - file (TFBS) from known_TFBS_for_overlap + set file (bed_footprints), file (bed_motifs), file (fasta) from for_overlap output: file ('*.FASTA') into FASTA_for_clustering script: """ + ${path_bin}/compareBed.sh --data ${bed_footprints} --motifs ${bed_motifs} --fasta ${fasta} -o ${name_placeholder} -min ${params.min_size_fp} -max ${params.max_size_fp} """ } @@ -53,7 +61,7 @@ process overlap_with_known_TFBS { process clustering { input: - file (faste) from FASTA_for_clustering + file (fasta) from FASTA_for_clustering output: set name, file ('*.bed') into bed_for_motif_esitmation @@ -78,7 +86,7 @@ process bed_to_clustered_fasta { script: """ - Rscript ${path_bin}/bed_to_fasta.R ${bed} ${name} + Rscript ${path_bin}/bed_to_fasta.R ${bed} ${name} ${params.min_seq} """ } @@ -101,7 +109,7 @@ process glam2 { script: """ - glam2 n ${fasta} -O . + glam2 n ${fasta} -O . -a 10 -b 20 -z 5 """ } @@ -122,12 +130,12 @@ process tomtom { script: """ - tomtom ${meme} ${jaspar_db} -thresh 0.01 -mi 1 -text | sed '/^#/ d' | sed '/^\$/d' > ${name}_known_motif.tsv + tomtom ${meme} ${jaspar_db} -thresh 0.01 -text --norc | sed '/^#/ d' | sed '/^\$/d' > ${name}_known_motif.tsv """ } -//Joining channels with meme and tsv files. Filter joined channel on linecount. +//Joining channels with meme and tsv files. Filter joined channel on line count. //Only meme-files which corresponding tsv files have linecount <= 1 are writen to next channel. for_filter = meme_for_filter.join( tsv_for_filter ) for_filter @@ -138,7 +146,7 @@ for_filter .into { meme_for_scan; check } -//If unknown channel is empty print errormessage +//If channel 'check' is empty print errormessage process check_for_unknown_motifs { echo true @@ -155,6 +163,7 @@ process check_for_unknown_motifs { } +//Get the best(first) Motif from each MEME-file process get_best_motif { input: @@ -212,15 +221,25 @@ process create_GTF { } -process create_UROPA_config { +bed_for_final_filter.combine(gtf_for_uropa).set {uropa_in} + + +// Create configuration file for UROPA. +// Takes template and replaces bed- and gtf-placeholders with actual paths. +process create_uropa_config { + + publishDir '/mnt/agnerds/Rene.Wiegandt/10_Master/', mode: 'copy' input: + set val(bed), val(gtf) from uropa_in.toList() + file (conf) from config output: file ('uropa.config') into uropa_config script: """ + sed -- 's/placeholder_gtf/${gtf}/g; s/placeholder_bed/${bed}/g' ${conf} > uropa.config.final """ } diff --git a/uropa.config b/uropa.config new file mode 100644 index 0000000..3c4a189 --- /dev/null +++ b/uropa.config @@ -0,0 +1,16 @@ +{ +"queries":[ + {"feature":"", + "feature.anchor": "", + "distance": [ , ], + "strand":"", + "direction":"", + "internals":"", + "filter.attribute":"", + "attribute.value":"", + "show.attributes":"" } + ], +"priority": "", +"gtf": "placeholder_gtf", +"bed": "placeholder_bed" +} \ No newline at end of file