Run gffread in multi-thread mode

Question

Is there any possibility to run gffread in multi-thread mode? The answer seems to be 'no' from the manual (or gffread -h), as no multi-thread option is mentioned.

I'm mostly using this utility to extract transcript sequences (FASTA) from annotation files (GTF). In single-thread mode, runtime is acceptable in most cases, but would like to speed up a bit in some cases (especially for big genomes).

Note that I already tried, as mentioned in the manual, to index the genome first (which, again, cuts the runtime in most but not all cases):

Note that the retrieval of the transcript sequences this way is going to be much faster if a fasta index file (genome.fa.fai in this example) is found in the same directory with the genomic fasta file.

Such an index file can be created with the samtools utility prior to running gffread, like this: samtools faidx genome.fa

Then in subsequent runs using the -g option gffread will find this FASTA index and use it to speed up the extraction of transcript sequences.

I'm also open to different alternatives than gffread to convert GTF -> FASTA.

Answer 1

Since there doesn't seem to be an easy way to run this in parallel, you could instead break the job into sections. For example, separate each chromosome into its own gtf file, extract the sequences using that file and then cat them all together.

The commands would be something like this (using the human hg38 genome and the gff GENCODE annotations downloaded from here):

Split the file into one file per chromosome (or whatever else you have as the first field of your gff/gtf file):

tmp=$(mktemp -d); ## create a temp directory
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf 

Extract the sequences described in each file, running each command in the background (&):

for f in $tmp/*; do 
    gffread -w "$f.fa" -g hg38.fa "$f" & 
done

Concatenate them into one file:

cat "$tmp/*fa" > all.fa

You can combine all this into a single command, using wait to make sure you don't concatenate until all of the sub-commands have finished:

tmp=$(mktemp -d); 
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf && 
    for f in "$tmp"/*; do 
        gffread-0.9.12.Linux_x86_64/gffread -w $f.fa -g hg38.fa "$f" & 
    done 
    wait 
    cat "$tmp"/*fa > all.fa

Or, just save it as a script called pgffread.sh (or whatever):

#!/bin/sh
gtf="$1"
genome="$2"

tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' "$gtf" &&
    for f in "$tmp"/*; do
      gffread -w "$f".fa -g "$genome" "$f" &
    done
wait
cat "$tmp"/*fa 
rm -rf "$tmp"

And then run it giving the genome and annotation files as parameters:

pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa

On the system I ran it on, which has 48 cores, the time difference was significant (although it's so fast anyway, I really don't know if it is worth the effort unless you do this very often):

$ time gffread -w out.fa -g hg38.fa gencode.v27.annotation.gtf
real    1m52.072s
user    1m3.160s
sys     0m47.257s

$ time ./pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa
real    0m17.533s
user    1m5.228s
sys     0m43.039s