Skip to content

Commit

Permalink
Give up on Jaccard indices
Browse files Browse the repository at this point in the history
Computations still crashed because the chromosome sort order should be
the same in the genomefile and in the BED file.  At this point I prefer
to just remove the whole computation.  Jaccard indices were not very
useful when one side is much wider than the other.
  • Loading branch information
charles-plessy committed Oct 9, 2024
1 parent 9159900 commit b46cfd4
Show file tree
Hide file tree
Showing 6 changed files with 0 additions and 70 deletions.
3 changes: 0 additions & 3 deletions modules/local/merge_repeatmasker_all.nf
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ process MERGE_REPM_RESULTS {

output:
tuple val(meta), path("*.fasta.gz") , emit: fasta
tuple val(meta), path("*_jaccard.txt") , emit: txt
tuple val(meta), path("*.bed.gz") , emit: bed_gz
tuple val(meta), path("*_repeatmasker_all.mask.bed.gz") , emit: repm_all_bed_gz
path "versions.yml" , emit: versions
Expand All @@ -24,9 +23,7 @@ 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" -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"
}
Expand Down
3 changes: 0 additions & 3 deletions modules/local/merge_repeatmasker_dfam.nf
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ process MERGE_REPM_RESULTS {

output:
tuple val(meta), path("*.fasta.gz") , emit: fasta
tuple val(meta), path("*_jaccard.txt") , emit: txt
tuple val(meta), path("*.bed.gz") , emit: bed_gz
tuple val(meta), path("*_repeatmasker_all.mask.bed.gz") , emit: repm_all_bed_gz
path "versions.yml" , emit: versions
Expand All @@ -24,9 +23,7 @@ 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" -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"
}
Expand Down
3 changes: 0 additions & 3 deletions modules/local/merge_repeatmasker_extlib.nf
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ process MERGE_REPM_RESULTS {

output:
tuple val(meta), path("*.fasta.gz") , emit: fasta
tuple val(meta), path("*_jaccard.txt") , emit: txt
tuple val(meta), path("*.bed.gz") , emit: bed_gz
tuple val(meta), path("*_repeatmasker_all.mask.bed.gz") , emit: repm_all_bed_gz
path "versions.yml" , emit: versions
Expand All @@ -24,9 +23,7 @@ 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" -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"
}
Expand Down
6 changes: 0 additions & 6 deletions modules/local/mergemasks.nf
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ process MERGE_MASKS {

output:
tuple val(meta), path("*.fasta.gz") , emit: fasta
tuple val(meta), path("*_jaccard.txt") , emit: txt
tuple val(meta), path("*.bed.gz") , emit: bed_gz
path "versions.yml" , emit: versions

Expand All @@ -23,11 +22,6 @@ process MERGE_MASKS {
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!
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
zcat $windowmasker $repeatmasker | sort -k1,1 -k2,2n | bedtools merge | gzip --best --no-name > ${prefix}_windowmasker_repeatmasker.bed.gz
Expand Down
48 changes: 0 additions & 48 deletions modules/local/multiqc_softmask_overlaps.nf

This file was deleted.

7 changes: 0 additions & 7 deletions workflows/pairgenomealignmask.nf
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,6 @@ include { MERGE_MASKS as MERGEDMASKS_ALL } from '../modules/
include { GFASTATS as MERGEDMASKS_ALL_STATS } from '../modules/nf-core/gfastats/main'

include { MULTIQC_SOFTMASK_STATS } from '../modules/local/multiqc_softmask_statistics.nf'
include { MULTIQC_SOFTMASK_OVERLAPS } from '../modules/local/multiqc_softmask_overlaps.nf'

include { MULTIQC } from '../modules/nf-core/multiqc/main'
include { paramsSummaryMap } from 'plugin/nf-validation'
Expand Down Expand Up @@ -170,12 +169,6 @@ workflow PAIRGENOMEALIGNMASK {
)
ch_multiqc_files = ch_multiqc_files.mix(MULTIQC_SOFTMASK_STATS.out.tsv)

// Aggregation of statistics (Jaccard indices)
//
MULTIQC_SOFTMASK_OVERLAPS ( MERGEDMASKS_ALL.out.txt.map{it[1]}.collect() )
ch_multiqc_files = ch_multiqc_files.mix(MULTIQC_SOFTMASK_OVERLAPS.out.tsv)


// Collect software versions
//
ch_versions = ch_versions
Expand Down

0 comments on commit b46cfd4

Please sign in to comment.