6
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ I have used gffread a few times. My impression is it is fast. As I have just tried, for Ensembl human annotations, it gets all transcripts in 29 sec. This should be faster than many other operations on human genome. How often do you want to run gffread? $\endgroup$
    – user172818
    Mar 14 '18 at 16:43
  • $\begingroup$ Thank you. Yes it is usually quite fast. Actually I was adressing this question since it took me more than 8min to get the transcripts out of a putative genome size of 1.4Gb, but could not reproduce the slow runtime on another server (where it took me around 110sec). Weird. $\endgroup$
    – aechchiki
    Mar 14 '18 at 21:40
3
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ thanks! yes, I should have thought of breaking the reference into smaller bits. $\endgroup$
    – aechchiki
    Mar 15 '18 at 7:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.