SEQanswers

Go Back   SEQanswers > Sequencing Technologies/Companies > Illumina/Solexa



Similar Threads
Thread Thread Starter Forum Replies Last Post
A question about BWA index louis7781x Bioinformatics 9 04-02-2012 03:23 AM
BWA error and reference index satishkumar Introductions 1 11-19-2010 08:22 AM
bwa index option bwtsw sowmyai Bioinformatics 2 07-15-2010 10:30 AM
BWA index error aldrinyim Bioinformatics 3 02-14-2010 11:03 AM
BWA index adrian Bioinformatics 3 01-23-2010 08:10 AM

Reply
 
Thread Tools
Old 12-03-2013, 04:28 AM   #21
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default bwa index

BWA users- I get less than 8 output files when using BWA to index the reference genome- any explanation, please and would this affect BWA alignment?

bwa index reference.fna
mmmm is offline   Reply With Quote
Old 12-03-2013, 04:53 AM   #22
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

Quote:
Originally Posted by mmmm View Post
BWA users- I get less than 8 output files when using BWA to index the reference genome- any explanation, please and would this affect BWA alignment?

bwa index reference.fna
Please post additional information about what OS (32/64-bit) you are using (virtual or real), amount of memory you have. Are there any error messages? I assume "reference.fna" is a multi-fasta file? Which files are you seeing (extensions)?
GenoMax is offline   Reply With Quote
Old 12-03-2013, 05:00 AM   #23
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

the reference.fna is not a multi-fasta files. it is the complete bacterial genome.
I can see (.amb, .ann, .bwt, .pac, .sa) but can not see (.rpac, .rsa, .rbwt)

I have more than enough space on the server. there is no error message?
mmmm is offline   Reply With Quote
Old 12-03-2013, 05:20 AM   #24
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

That looks right.

As of bwa version 0.6.x the forward and reverse indexes were merged, eliminating the .rxxx files.
GenoMax is offline   Reply With Quote
Old 12-03-2013, 05:50 AM   #25
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

grateful thanks for clarification. Now, I am trying to detect the gene coverage cin the mapped reads compared to the reference genome using

bedtools
coverageBed -abam alignment.bam -b gff > coverage.txt

(but I have found that no genes have been mapped at all to the reference genome when opened the coverage.txt in a libreoffice)?
mmmm is offline   Reply With Quote
Old 12-03-2013, 07:20 AM   #26
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

See the notes for "genomecov": http://bedtools.readthedocs.org/en/l...genomecov.html

Have you sorted your bam file?
GenoMax is offline   Reply With Quote
Old 12-03-2013, 09:52 AM   #27
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

I did sort the bam file- should this be ok?

bedtools
coverageBed -abam sortedbam -b gff > coverage.txt
mmmm is offline   Reply With Quote
Old 12-03-2013, 12:09 PM   #28
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

Quote:
Originally Posted by mmmm View Post
I did sort the bam file- should this be ok?

bedtools
coverageBed -abam sortedbam -b gff > coverage.txt
That should be ok. Are you not seeing anything in the coverage.txt file?

Do you want to try "multicov" along with a BED format file for the genes? http://bedtools.readthedocs.org/en/l.../multicov.html
GenoMax is offline   Reply With Quote
Old 12-04-2013, 01:43 AM   #29
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

the coverage.txt file contains genes annotations (gff) but none of my aligned reads are covered by these annotated genes; I get (0) in the last column of the coverage.txt file for all genes
mmmm is offline   Reply With Quote
Old 12-04-2013, 04:03 AM   #30
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

The "gff" in your command line above is the actual GFF format file for the annotations, correct (here if the format reference for GFF3: http://www.sequenceontology.org/gff3.shtml And the file actually matches the original genome used for alignment?

Have you checked your aligned bam file (using IGV for example) to make sure that the alignments look good?

Last edited by GenoMax; 12-04-2013 at 04:24 AM.
GenoMax is offline   Reply With Quote
Old 12-04-2013, 04:13 AM   #31
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

yes, the alignment looks very good using IGV- I think, there is an issue in the GFF file as genes do not appear annotated using IGV (only the sequence appear without gene annotation), do not know how to fix the gff file?
mmmm is offline   Reply With Quote
Old 12-04-2013, 04:14 AM   #32
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

Can you post an example of the GFF file (first few lines would be fine)?
GenoMax is offline   Reply With Quote
Old 12-04-2013, 04:22 AM   #33
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

>complete genome
NC_022544.1 RefSeq region 1 4814801 . + . ID=id0;Dbxref=taxon:568709;Is_circular=true;gbkey=Src;genome=genomic;mol_type=genomic DNA;serovar=Typhimurium;strain=DT2;sub-species=enterica
NC_022544.1 RefSeq gene 169 255 . + . ID=gene0;Name=thrL;Dbxref=GeneID:17155329;gbkey=Gene;gene=thrL;locus_tag=STMDT2_00011
NC_022544.1 RefSeq CDS 169 255 . + 0 ID=cds0;Name=YP_008642919.1;Parent=gene0;Dbxref=Genbank:YP_008642919.1,GeneID:17155329;gbkey=CDS;product=thr operon leader peptide;protein_id=YP_008642919.1;transl_table=11
NC_022544.1 RefSeq gene 337 2799 . + . ID=gene1;Name=thrA;Dbxref=GeneID:17159252;gbkey=Gene;gene=thrA;locus_tag=STMDT2_00021
NC_022544.1 RefSeq CDS 337 2799 . + 0 ID=cds1;Name=YP_008642920.1;Parent=gene1;Dbxref=Genbank:YP_008642920.1,GeneID:17159252;gbkey=CDS;product=aspartokinase I%2Fhomoserine dehydrogenase I;protein_id=YP_008642920.1;transl_table=11
NC_022544.1 RefSeq gene 2801 3730 . + . ID=gene2;Name=thrB;Dbxref=GeneID:17159249;gbkey=Gene;gene=thrB;locus_tag=STMDT2_00031
NC_022544.1 RefSeq CDS 2801 3730 . + 0 ID=cds2;Name=YP_008642921.1;Parent=gene2;Dbxref=Genbank:YP_008642921.1,GeneID:17159249;gbkey=CDS;product=Homoserine kinase;protein_id=YP_008642921.1;transl_table=11
NC_022544.1 RefSeq gene 3734 5020 . + . ID=gene3;Name=thrC;Dbxref=GeneID:17159250;gbkey=Gene;gene=thrC;locus_tag=STMDT2_00041
NC_022544.1 RefSeq CDS 3734 5020 . + 0 ID=cds3;Name=YP_008642922.1;Parent=gene3;Dbxref=Genbank:YP_008642922.1,GeneID:17159250;gbkey=CDS;product=threonine synthase;protein_id=YP_008642922.1;transl_table=11
NC_022544.1 RefSeq gene 5114 5887 . - . ID=gene4;Name=yaaA;Dbxref=GeneID:17159251;gbkey=Gene;gene=yaaA;locus_tag=STMDT2_00051
NC_022544.1 RefSeq CDS 5114 5887 . - 0 ID=cds4;Name=YP_008642923.1;Parent=gene4;Dbxref=Genbank:YP_008642923.1,GeneID:17159251;gbkey=CDS;product=hypothetical protein;protein_id=YP_008642923.1;transl_table=11
NC_022544.1 RefSeq gene 5966 7396 . - . ID=gene5;Name=yaaJ;Dbxref=GeneID:17159391;gbkey=Gene;gene=yaaJ;locus_tag=STMDT2_00061
NC_022544.1 RefSeq CDS 5966 7396 . - 0 ID=cds5;Name=YP_008642924.1;Parent=gene5;Dbxref=Genbank:YP_008642924.1,GeneID:17159391;gbkey=CDS;product=putative amino-acid transport protein;protein_id=YP_008642924.1;transl_table=11
NC_022544.1 RefSeq gene 7665 8618 . + . ID=gene6;Name=talB;Dbxref=GeneID:17159395;gbkey=Gene;gene=talB;locus_tag=STMDT2_00071
NC_022544.1 RefSeq CDS 7665 8618 . + 0 ID=cds6;Name=YP_008642925.1;Parent=gene6;Dbxref=Genbank:YP_008642925.1,GeneID:17159395;gbkey=CDS;product=transaldolase B;protein_id=YP_008642925.1;transl_table=11
NC_022544.1 RefSeq gene 8729 9319 . + . ID=gene7;Name=mog;Dbxref=GeneID:17159215;gbkey=Gene;gene=mog;locus_tag=STMDT2_00081
NC_022544.1 RefSeq CDS 8729 9319 . + 0 ID=cds7;Name=YP_008642926.1;Parent=gene7;Dbxref=Genbank:YP_008642926.1,GeneID:17159215;gbkey=CDS;product=molybdopterin biosynthesis Mog protein;protein_id=YP_008642926.1;transl_table=11
NC_022544.1 RefSeq gene 9376 9942 . - . ID=gene8;Name=yaaH;Dbxref=GeneID:17158379;gbkey=Gene;gene=yaaH;locus_tag=STMDT2_00091
NC_022544.1 RefSeq CDS 9376 9942 . - 0
mmmm is offline   Reply With Quote
Old 12-04-2013, 04:34 AM   #34
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

Can you remove the fist line from the file (please make a backup copy of original file in case something goes wrong)

Code:
>complete genome
and then try? That would make it a gff format file (http://www.sanger.ac.uk/resources/so.../gff/spec.html).

To make it a gff3 format file you will have to replace that first line with following two lines http://www.sequenceontology.org/gff3.shtml

Code:
##gff-version 3 
##sequence-region NC_022544.1 1 4814801

Last edited by GenoMax; 12-04-2013 at 06:08 AM.
GenoMax is offline   Reply With Quote
Old 12-04-2013, 05:08 AM   #35
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

after editing the first 2 lines in the gff file as you have kindly suggested. genes annotations can not be seen on IGV (can see only the sequence but not the genes)???- your advice is very appreciated
mmmm is offline   Reply With Quote
Old 12-04-2013, 05:26 AM   #36
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

Have you renamed the GFF file as "your_file_name.gff3"?

IGV expects the GFF3 files to have that extension (http://www.broadinstitute.org/software/igv/GFF) and they also need to be tab-delimited (which your example above does not appear to be).
GenoMax is offline   Reply With Quote
Old 12-04-2013, 06:28 AM   #37
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

After taking out the first two lines from your example, adding the two lines for GFF3 meta-data and then converting to tab-delimited text the file appears to work with IGV.

Take these two lines out:

Code:
>complete genome
NC_022544.1 RefSeq region 1 4814801 . + . ID=id0;Dbxref=taxon:568709;Is_circular=true;gbkey=Src;genome=genomic;mol_type=genomic DNA;serovar=Typhimurium;strain=DT2;sub-species=enterica
Replace with:

Code:
##gff-version 3 
##sequence-region NC_022544.1 1 4814801
Attached Images
File Type: png IGV_Cap.PNG (14.2 KB, 2 views)
GenoMax is offline   Reply With Quote
Old 12-05-2013, 02:11 AM   #38
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

how do you convert text (gff3) to tab-delimited, please?
I used:

awk '{ for(i=1;i<=NF;i++){if(i==NF){printf("%s\n",$NF);}else {printf("%s\t",$i)}}}' file.gff3

but did not convert the file to tab-delimited????

##gff-version 3
##sequence-region NC_022544.1 1 4814801

NC_022544.1 RefSeq gene 169 255 . + . ID=gene0;Name=thrL;Dbxref=GeneID:17155329;gbkey=Gene;gene=thrL;locus_tag=STMDT2_00011
NC_022544.1 RefSeq CDS 169 255 . + 0 ID=cds0;Name=YP_008642919.1;Parent=gene0;Dbxref=Genbank:YP_008642919.1,GeneID:17155329;gbkey=CDS;gene=thrL;product=thr operon leader peptide;protein_id=YP_008642919.1;transl_table=11
NC_022544.1 RefSeq gene 337 2799 . + . ID=gene1;Name=thrA;Dbxref=GeneID:17159252;gbkey=Gene;gene=thrA;locus_tag=STMDT2_00021
NC_022544.1 RefSeq CDS 337 2799 . + 0 ID=cds1;Name=YP_008642920.1;Parent=gene1;Dbxref=Genbank:YP_008642920.1,GeneID:17159252;gbkey=CDS;gene=thrA;product=aspartokinase I%2Fhomoserine dehydrogenase I;protein_id=YP_008642920.1;transl_table=11
NC_022544.1 RefSeq gene 2801 3730 . + . ID=gene2;Name=thrB;Dbxref=GeneID:17159249;gbkey=Gene;gene=thrB;locus_tag=STMDT2_00031
NC_022544.1 RefSeq CDS 2801 3730 . + 0 ID=cds2;Name=YP_008642921.1;Parent=gene2;Dbxref=Genbank:YP_008642921.1,GeneID:17159249;gbkey=CDS;gene=thrB;product=Homoserine kinase;protein_id=YP_008642921.1;transl_table=11
NC_022544.1 RefSeq gene 3734 5020 . + . ID=gene3;Name=thrC;Dbxref=GeneID:17159250;gbkey=Gene;gene=thrC;locus_tag=STMDT2_00041
NC_022544.1 RefSeq CDS 3734 5020 . + 0 ID=cds3;Name=YP_008642922.1;Parent=gene3;
mmmm is offline   Reply With Quote
Old 12-05-2013, 04:36 AM   #39
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,075
Default

Try this unix command. Adjust the file names accordingly.
Code:
$ tr ' ' \\t < original.gff3 > tab_converted.gff3
There is a "space" between the two single quotes in the command above.

Best to put the two top metadata lines in after the conversion.
GenoMax is offline   Reply With Quote
Old 12-05-2013, 05:35 AM   #40
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default

thanks you soo much for your advice and time. managed to view genes on IGV and have got a coverage.txt file showing genes that are covered/ missed from the reference

but have a very simple technichal issue when I open the coverage.txt using excel or libreoffice- I could not sort the coverage values in ascending order as values are written on a separte line

NC_022544.1 RefSeq gene 2096621 2097676 . - . ID=gene2037;Name=cbiG;Dbxref=GeneID:17157414;gbkey=Gene;gene=cbiG;locus_tag=STMDT2_20011
2505 1056 1056 1
NC_022544.1 RefSeq CDS 2096621 2097676 . - 0 "ID=cds1963;Name=YP_008644883.1;Parent=gene2037;Dbxref=Genbank:YP_008644883.1,GeneID:17157414;gbkey=CDS;gene=cbiG;product=cobalamin biosynthesis protein;protein_id=YP_008644883.1;transl_table=11"
2505 1056 1056 0.9
NC_022544.1 RefSeq gene 3145124 3146200 . - . ID=gene2977;Name=STMDT2_29251;Dbxref=GeneID:17156485;gbkey=Gene;locus_tag=STMDT2_29251
1716 1077 1077 0.6
mmmm is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:13 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO