RSEM implementation

Question

[anasjamshed]

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?

[rpolicastro]

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.

[GenoMax]

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

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

[anasjamshed]

which index?

[Ram]

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.

[anasjamshed]

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!

[GenoMax]

Error is clearly noted. Check your gff file :

Line 25 does not have 9 fields

[anasjamshed]

It contains 24 lines only and does not have 25 line

[Ram]

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.

[anasjamshed]

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?

Answer 1

[Ram]

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

[anasjamshed]

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!

[Ram]

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.

[anasjamshed]

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

[Ram]

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.