This is a test version of Biostars. For the public version, visit https://www.biostars.org.
RSEM implementation

I have the virus genome(fasta) and gff file and I am trying to prepare-reference through the following commands:

rsem-prepare-reference --gff3 KT992094.1.gff3 KT992094.1.fasta

or

rsem-prepare-reference --gff3 KT992094.1.gff \ --gff3-genes-as-transcripts \ --bowtie \ KT992094.1.fasta \ ref/virus

But it's saying:

Invalid number of arguments!

How can I solve this issue?

rsem
• 4,706 views • link
updated by Ram Ram • written by anasjamshed anasjamshed

Try removing the backslashes and extra spaces in your second command. Backslashes escape characters, so if it's all one one line you'll be escaping important characters in your arguments and flags.

Also make sure that you're using the correct command line arguments and flags, since your two commands vary greatly in what they theoretically are doing.

• 0 views • link
written by rpolicastro rpolicastro 
rsem-prepare-reference --gff3 KT992094.1.gff3 KT992094.1.fasta basename

basesname will be the name for resulting index etc. So provide that.

• 0 views • link
written by GenoMax GenoMax

The index you're trying to generate. Ensure you're either in the directory you wish the output of rsem-prepare-reference to be in, or you give that output directory along with a sensible prefix as the basename.

• 0 views • link
written by Ram Ram

while running this:

rsem-prepare-reference --gff3 KT992094.1.gff3 KT992094.1.fasta out

I got this error:

rsem-gff3-to-gtf KT992094.1.gff3 out.gtf
Line 25 does not have 9 fields:

"rsem-gff3-to-gtf KT992094.1.gff3 out.gtf" failed! Plase check if you provide correct parameters/options for the pipeline!
• 0 views • link
written by anasjamshed anasjamshed

Error is clearly noted. Check your gff file :

Line 25 does not have 9 fields

• 0 views • link
written by GenoMax GenoMax

gff3 file

It contains 24 lines only and does not have 25 line

• 0 views • link
written by anasjamshed anasjamshed

Are you sure? Paste the output to cat -te KT992094.1.gff3 && echo $? here - don't take a screenshot, copy and paste the content and code-format it.

• 0 views • link
written by Ram Ram

I tried :

cat -te KT992094.1.gff3 && echo $?

And get the following results:

##sequence-region KT992094.1 1 15223$
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=11250$
KT992094.1^IGenbank^Iregion^I1^I15223^I.^I+^I.^IID=KT992094.1:1..15223;Dbxref=taxon:11250;gb-acronym=HRSV;gbkey=Src;genome=genomic;mol_type=viral cRNA;note=recombinant D46/D53 strain;old-name=Human respiratory syncytial virus;strain=A2$
KT992094.1^IGenbank^Igene^I45^I576^I.^I+^I.^IID=gene-NS1;Name=NS1;gbkey=Gene;gene=NS1;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I99^I518^I.^I+^I0^IID=cds-ALS35583.1;Parent=gene-NS1;Dbxref=NCBI_GP:ALS35583.1;Name=ALS35583.1;gbkey=CDS;gene=NS1;product=nonstructural protein 1;protein_id=ALS35583.1$
KT992094.1^IGenbank^Igene^I596^I1098^I.^I+^I.^IID=gene-NS2;Name=NS2;gbkey=Gene;gene=NS2;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I628^I1002^I.^I+^I0^IID=cds-ALS35584.1;Parent=gene-NS2;Dbxref=NCBI_GP:ALS35584.1;Name=ALS35584.1;gbkey=CDS;gene=NS2;product=nonstructural protein 2;protein_id=ALS35584.1$
KT992094.1^IGenbank^Igene^I1126^I2328^I.^I+^I.^IID=gene-N;Name=N;gbkey=Gene;gene=N;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I1141^I2316^I.^I+^I0^IID=cds-ALS35585.1;Parent=gene-N;Dbxref=NCBI_GP:ALS35585.1;Name=ALS35585.1;gbkey=CDS;gene=N;product=nucleoprotein;protein_id=ALS35585.1$
KT992094.1^IGenbank^Igene^I2330^I3243^I.^I+^I.^IID=gene-P;Name=P;gbkey=Gene;gene=P;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I2347^I3072^I.^I+^I0^IID=cds-ALS35586.1;Parent=gene-P;Dbxref=NCBI_GP:ALS35586.1;Name=ALS35586.1;gbkey=CDS;gene=P;product=phosphoprotein;protein_id=ALS35586.1$
KT992094.1^IGenbank^Igene^I3253^I4210^I.^I+^I.^IID=gene-M;Name=M;gbkey=Gene;gene=M;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I3262^I4032^I.^I+^I0^IID=cds-ALS35587.1;Parent=gene-M;Dbxref=NCBI_GP:ALS35587.1;Name=ALS35587.1;gbkey=CDS;gene=M;product=matrix protein;protein_id=ALS35587.1$
KT992094.1^IGenbank^Igene^I4220^I4629^I.^I+^I.^IID=gene-SH;Name=SH;gbkey=Gene;gene=SH;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I4304^I4498^I.^I+^I0^IID=cds-ALS35582.1;Parent=gene-SH;Dbxref=NCBI_GP:ALS35582.1;Name=ALS35582.1;gbkey=CDS;gene=SH;product=small hydrophobic protein;protein_id=ALS35582.1$
KT992094.1^IGenbank^Igene^I4674^I5596^I.^I+^I.^IID=gene-G;Name=G;gbkey=Gene;gene=G;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I4689^I5585^I.^I+^I0^IID=cds-ALS35588.1;Parent=gene-G;Dbxref=NCBI_GP:ALS35588.1;Name=ALS35588.1;gbkey=CDS;gene=G;product=attachment glycoprotein G;protein_id=ALS35588.1$
KT992094.1^IGenbank^Igene^I5649^I7551^I.^I+^I.^IID=gene-F;Name=F;gbkey=Gene;gene=F;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I5662^I7386^I.^I+^I0^IID=cds-ALS35589.1;Parent=gene-F;Dbxref=NCBI_GP:ALS35589.1;Name=ALS35589.1;gbkey=CDS;gene=F;product=fusion protein;protein_id=ALS35589.1$
KT992094.1^IGenbank^Igene^I7598^I8558^I.^I+^I.^IID=gene-M2;Name=M2;gbkey=Gene;gene=M2;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I7607^I8191^I.^I+^I0^IID=cds-ALS35591.1;Parent=gene-M2;Dbxref=NCBI_GP:ALS35591.1;Name=ALS35591.1;gbkey=CDS;gene=M2;product=m2-1;protein_id=ALS35591.1$
KT992094.1^IGenbank^ICDS^I8160^I8432^I.^I+^I0^IID=cds-ALS35592.1;Parent=gene-M2;Dbxref=NCBI_GP:ALS35592.1;Name=ALS35592.1;gbkey=CDS;gene=M2;product=m2-2 protein;protein_id=ALS35592.1$
KT992094.1^IGenbank^Igene^I8491^I15068^I.^I+^I.^IID=gene-L;Name=L;gbkey=Gene;gene=L;gene_biotype=protein_coding$
KT992094.1^IGenbank^ICDS^I8499^I14996^I.^I+^I0^IID=cds-ALS35590.1;Parent=gene-L;Dbxref=NCBI_GP:ALS35590.1;Name=ALS35590.1;gbkey=CDS;gene=L;product=L polymerase protein;protein_id=ALS35590.1$
$
0

This file does not contain 25th line?

• 0 views • link
written by anasjamshed anasjamshed

1 answer

Your GFF3 file has 25th line - see the lone $ in the cat -te output? That's an empty line. Delete the empty line using

sed '/^$/d' KT992094.1.gff3 > KT992094.1.new.gff3

And then run your old command (with the new gff3 file)

rsem-prepare-reference --gff3 KT992094.1.new.gff3 KT992094.1.fasta out
• 0 views • link
written by Ram Ram

Thanks. The command ran successfully but it threw the following errors:

rsem-gff3-to-gtf KT992094.1.new.gff3 out.gtf
GTF file is successully generated.
There are 0 transcripts contained in the generated GTF file.

rsem-extract-reference-transcripts out 0 out.gtf None 0 KT992094.1.fasta
The reference contains no transcripts!
"rsem-extract-reference-transcripts out 0 out.gtf None 0 KT992094.1.fasta" failed! Plase check if you provide correct parameters/options for the pipeline!
• 0 views • link
written by anasjamshed anasjamshed

The error is pretty self-descriptive. Where did you get this gff3 file? You need a better transcript feature file. Also, this is a new topic so you should open a new question.

I'll move my comment to an answer, please accept it to provide closure to this post.

• 0 views • link
written by Ram Ram

I run the following command, and it successfully finds the transcripts:

 rsem-prepare-reference --gff3 KT992094.1.new.gff3 --gff3-genes-as-transcripts --bowtie KT992094.1.fasta rsv-rsem
• 0 views • link
written by anasjamshed anasjamshed

That's wonderful - congratulations on finding the solution yourself, and thank you for sharing it here. Please accept my answer to provide closure to the post.

• 0 views • link
written by Ram Ram

Log in to answer this question.

More posts like this