Create a "genome file" that bedtools jaccard needs.

Otherwise it crashes if a sequence from B is not in A. There is obvious repetition of the same computation here. It is not satisfying but I need a short term fix. Hopefully it will put the genome sequence in the node's memory cache, which will be useful in the final commands.
oist · Oct 9, 2024 · 9159900 · 9159900
1 parent 264035e
commit 9159900
Show file tree

Hide file tree

Showing 4 changed files with 10 additions and 6 deletions.
diff --git a/modules/local/merge_repeatmasker_all.nf b/modules/local/merge_repeatmasker_all.nf
@@ -24,8 +24,9 @@ process MERGE_REPM_RESULTS {
     def args = task.ext.args ?: ''
     def prefix = task.ext.prefix ?: "${meta.id}"
     """
+    awk '/^>/ {if (seqlen){print seqname "\t" seqlen}; split(\$1, a, ">"); seqname=a[2]; seqlen=0; next} {seqlen += length(\$0)} END {print seqname "\t" seqlen}' $genome > genome.genome # thanks, ChatGPT!
     run_bedtools_operations() {
-        bedtools jaccard -nonamecheck -a "\$1" -b "\$2"                               > "${prefix}_\${3}_jaccard.txt"
+        bedtools jaccard -nonamecheck -a "\$1" -b "\$2" -g genome.genome              > "${prefix}_\${3}_jaccard.txt"
         zcat "\$1" "\$2" | sort -k1,1 -k2,2n | bedtools merge | gzip --best --no-name > "${prefix}_\${3}.mask.bed.gz"
     }
 

diff --git a/modules/local/merge_repeatmasker_dfam.nf b/modules/local/merge_repeatmasker_dfam.nf
@@ -24,8 +24,9 @@ process MERGE_REPM_RESULTS {
     def args = task.ext.args ?: ''
     def prefix = task.ext.prefix ?: "${meta.id}"
     """
+    awk '/^>/ {if (seqlen){print seqname "\t" seqlen}; split(\$1, a, ">"); seqname=a[2]; seqlen=0; next} {seqlen += length(\$0)} END {print seqname "\t" seqlen}' $genome > genome.genome # thanks, ChatGPT!
     run_bedtools_operations() {
-        bedtools jaccard -nonamecheck -a "\$1" -b "\$2"                               > "${prefix}_\${3}_jaccard.txt"
+        bedtools jaccard -nonamecheck -a "\$1" -b "\$2" -g genome.genome              > "${prefix}_\${3}_jaccard.txt"
         zcat "\$1" "\$2" | sort -k1,1 -k2,2n | bedtools merge | gzip --best --no-name > "${prefix}_\${3}.mask.bed.gz"
     }
 

diff --git a/modules/local/merge_repeatmasker_extlib.nf b/modules/local/merge_repeatmasker_extlib.nf
@@ -24,8 +24,9 @@ process MERGE_REPM_RESULTS {
     def args = task.ext.args ?: ''
     def prefix = task.ext.prefix ?: "${meta.id}"
     """
+    awk '/^>/ {if (seqlen){print seqname "\t" seqlen}; split(\$1, a, ">"); seqname=a[2]; seqlen=0; next} {seqlen += length(\$0)} END {print seqname "\t" seqlen}' $genome > genome.genome # thanks, ChatGPT!
     run_bedtools_operations() {
-        bedtools jaccard -nonamecheck -a "\$1" -b "\$2"                               > "${prefix}_\${3}_jaccard.txt"
+        bedtools jaccard -nonamecheck -a "\$1" -b "\$2" -g genome.genome              > "${prefix}_\${3}_jaccard.txt"
         zcat "\$1" "\$2" | sort -k1,1 -k2,2n | bedtools merge | gzip --best --no-name > "${prefix}_\${3}.mask.bed.gz"
     }
 

diff --git a/modules/local/mergemasks.nf b/modules/local/mergemasks.nf
@@ -23,9 +23,10 @@ process MERGE_MASKS {
     def args = task.ext.args ?: ''
     def prefix = task.ext.prefix ?: "${meta.id}"
     """
-    bedtools jaccard -nonamecheck -a $tantan       -b $windowmasker         > ${prefix}_tantan_windowmasker_jaccard.txt
-    bedtools jaccard -nonamecheck -a $tantan       -b $repeatmasker         > ${prefix}_tantan_repeatmasker_jaccard.txt
-    bedtools jaccard -nonamecheck -a $repeatmasker -b $windowmasker         > ${prefix}_repeatmasker_windowmasker_jaccard.txt
+    awk '/^>/ {if (seqlen){print seqname "\t" seqlen}; split(\$1, a, ">"); seqname=a[2]; seqlen=0; next} {seqlen += length(\$0)} END {print seqname "\t" seqlen}' $genome > genome.genome # thanks, ChatGPT!
+    bedtools jaccard -nonamecheck -a $tantan       -b $windowmasker -g genome.genome > ${prefix}_tantan_windowmasker_jaccard.txt
+    bedtools jaccard -nonamecheck -a $tantan       -b $repeatmasker -g genome.genome > ${prefix}_tantan_repeatmasker_jaccard.txt
+    bedtools jaccard -nonamecheck -a $repeatmasker -b $windowmasker -g genome.genome > ${prefix}_repeatmasker_windowmasker_jaccard.txt
 
     zcat $tantan $windowmasker               | sort -k1,1 -k2,2n | bedtools merge | gzip --best --no-name > ${prefix}_tantan_windowmasker.bed.gz
     zcat $tantan $repeatmasker               | sort -k1,1 -k2,2n | bedtools merge | gzip --best --no-name > ${prefix}_tantan_repeatmasker.bed.gz