Hello! I'm trying to build a phylogeny. I have a main species that I'm working with that has a genome assembly and GFF file (Species 1, white-throated sparrow, Z. albicolis) and several closely related species with whole genome sequencing, but no assembly (for simplicity, Species 2, Harris sparrow, Z. querula). Using the GATK pipeline, I aligned the reads from Species 2 to the Species 1 reference genome, called SNPs, generated a VCF file and then I used the GATKAlternateReferenceMaker to create a reference genome with the SNPs replaced. I'm then trying to use agat_sp_extract_sequences.pl to extract coding sequences from this new reference genome, but this results in a file that includes the coding sequence name headers, but empty sequences, e.g.
>rna-XM_005493870.3 gene=gene-LOC102071259 name=LOC102071259 seq_id=NW_005081745.1 type=cds
>rna-XM_026797357.1 gene=gene-LOC102073079 name=LOC102073079 seq_id=NW_005081745.1 type=cds
>rna-XM_005493869.3 gene=gene-LRWD1 name=LRWD1 seq_id=NW_005081745.1 type=cds
Here is the code used to produce the alternate reference genome:
$ gatk FastaAlternateReferenceMaker \
-R GCF_000385455.1_Zonotrichia_albicollis-1.0.1_genomic.fna \
-O Harris_sparrow.fasta \
-V Zonotrichia_querula.snps.vcf \
The output produced here looks, as best as I can tell, normal and the size of Harris_sparrow.fasta is comparable to the original genome of Species 1 (1021MB for Species 2 vs 1018MB). Next I ran the agat_sp_extract_sequences.pl script using a .GFF file that includes only the longest coding sequences (I used this same .GFF file for the same thing using different individuals from Species 1 without issue).
$ conda activate agat1
$ agat_sp_extract_sequences.pl -g GCF_000385455.1_Zonotrichia_albicollis-1.0.1_genomic.longest.gff -f Harris_sparrow.fasta -t cds -o Harris_sparrow_genomic.cds.fa
This code runs, albeit with some errors. Most notably, I have many errors listing "WARNING: Problem ! ID NW_005081924.1 not found !", where these entries starting with NW_are scaffold names in the reference genome of Species 1. The code ends in "14385 cds converted in fasta.", but it produces a file with the empty coding sequences above.
I am unable to figure out what is causing the issue... I'm assuming there is some kind of formatting error in the files produced, but I haven't been able to figure it out. Any ideas??
•
written
by
nmbaran •
1 answer
Ok, for completeness, I'm posting the solution here:
#Removed line numbers so that the line started with >NW_... for each scaffold/contig
sed 's/>[^N]*/>/g' Harris_sparrow.fasta > Harris_sparrow_fixed.fasta
#Removed : and everything after
sed 's/:.*//' Harris_sparrow_fixed.fasta > Harris_sparrow_fixed2.fasta
#Normalized line numbers in the fixed fasta file
picard NormalizeFasta \
I=Harris_sparrow_fixed2.fasta \
O=normalized_Harris_sparrow_fixed2.fasta
#Checked that it worked
samtools faidx normalized_Harris_sparrow_fixed2.fasta
#Ran!
agat_sp_extract_sequences.pl -g ../GCF_000385455.1_Zonotrichia_albicollis-1.0.1_genomic.longest.gff -f normalized_Harris_sparrow_fixed2.fasta -t cds -o Harris_sparrow_genomic.cds.fa
Log in to answer this question.
-
Need help with AUGUSTUS and AGAT for gene prediction and feature extraction, respectively.
written by Vijith •
-
Why are gene models from JGI non-sense when using AGAT to to extract proteins?
written by O.rka
-
Splitting multifasta by string match
written by Winston •
-
How to obtain protein-coding sequences from assembled genome/exome dataset?
written by DNAngel
-
SnpEff START/STOP codon not found error
written by thomasbersez
-
Error wiht MUMMER mapview
written by DanielC
-
Get CDS from consensus genome assembly
written by YocelynGG •
-
How to extract whole gene sequence with 5' and 3' UTR from genome
written by Bioinfonext
-
Software for Mcdonald-Kreitman test
written by zhangtao •
-
How to get gene coding sequences based on GFF file?
written by biolab
Is it just that there is a '>' before your fasta file for the '-f' input for
agat_sp_extract_sequences.pl?Otherwise, based on the warning messages, it would suggest that the names of the contigs have changed.
No, that isn't in the actual code. I added the '>' to try to help with formatting when writing this question.
Would part of the GATK pipeline change the contig names? I did need to use picard to add reading groups, e.g. https://gatk.broadinstitute.org/hc/en-us/articles/360035532352-Errors-about-read-group-RG-information
Not sure why it would have changed (although probably during the FastaAlternateReferenceMaker step) but you can just check the headers using
grep '>' Harris_sparrow.fasta. Perhaps just see if 'NW_005081924.1' is thereI did that and it did find that supposedly non-existent contig/scaffold! It returned this:
Sorry but what is the '389'? and usually the contig:1-180953 like scheme is when the regions '1-180953' of 'contig' have been extracted Is it likely that this addition of the coordinates to the end of the contig name is the reason?
You could try by removing everything in your contig names after the ':'
POST-EDIT: Ok after the formatting modification it is clearer. The headers have been modified and you need to trim back the name to the original accession name e.g. NW_005081924.1
I may have figured out the issue (but not 100% how to solve it yet)--just as you identified the potential problem, too. I attempted to extract the sequences instead using gffread and got a similar error, which led to more useful Googling. It looks like there is a name issue with the fasta file produced by FastaAlternateReferenceMaker, which is that it adds ">" and a line number to each entry of my fasta in Harris_sparrow.fasta, e.g.
See also:
I think I need to remove these line numbers from the fasta file and then retry, as neither gffread or agat properly parse those fasta headers. I need to step away, but I have a good direction to try now.
(Thanks so much for your help!)
To clarify In fasta header
“>389 NW_005081924.1:1-180953”
The sequence identifier is 389 all the rest after the first space is description. So in your GFF the sequence Identifier stored into the first column is NW_005081924.1 while the identifier in the fasta file is 389, so it cannot recognize / sync the information.
For future reference, you can format the
codepart of your post by using101010button (after selecting the text). I have done this for you in current post.>is used forquotingparts of text when responding. Like I did above.