SEQanswers - Bioinformatics

samtools pileup output for SNP detection

Mansequencer — Wed, 08 Sep 2010 15:56:15 GMT

Hi All,
I aligned 50 bp solexa reads onto a reference genome using bowtie and now I am using samtools pileup -vcf to call SNPs. This is the output it generates:
scaffold179 29989 T A 30 30 60 1 A B
scaffold179 29991 A G 4 4 60 1 G #
scaffold179 29992 A T 4 4 60 1 T #
scaffold179 29993 T A 4 4 60 1 N #
scaffold179 29994 T A 4 4 60 1 N #
scaffold179 29995 A C 4 4 60 1 C #
scaffold179 29997 T A 4 4 60 1 N #
scaffold179 29998 G C 4 4 60 1 C #
scaffold179 29999 T C 4 4 60 1 C #
scaffold179 30000 T A 4 4 60 1 N #
scaffold179 30001 T A 4 4 60 1 A #
scaffold179 30002 G A 4 4 60 1 A #
scaffold179 30003 A G 4 4 60 1 G #
scaffold179 30004 T C 4 4 60 1 C #
scaffold179 30005 T A 4 4 60 1 A #
scaffold179 30006 T A 4 4 60 1 N #
scaffold179 30007 G A 4 4 60 1 N #
scaffold179 30008 T A 4 4 60 1 N #
scaffold179 30009 G A 4 4 60 1 N #
scaffold179 30010 T A 4 4 60 1 N #
scaffold179 30011 T A 4 4 60 1 N #
scaffold179 30012 T A 4 4 60 1 N #
scaffold179 30013 T A 4 4 60 1 N #
scaffold179 30014 G A 4 4 60 1 N #

I understand that columns from left to right are:
Scaffold name:
coordinate:
reference base:
polymorphic base:

I am confused about the numbers after the 4th column. Please suggest as to what they correspond to.
Also when I try to look at the alignment in tview, it shows only Ns after the first 80 pb. My reference genome does contain some short stretches of Ns but is not all Ns so why does the samtools tview shows only Ns?
Any suggestion would be highly appreciated.
Thanks for all the help.
Thanks.

Lower case characters in FASTa reference sequence

foxyg — Wed, 08 Sep 2010 15:52:20 GMT

Hi, I have noticed in HG19, there are sometimes a segment of sequence that are in lower case? Does this have a special meaning?

Thanks

How to create a SAM/BAM file from scratch

DNAjunk — Wed, 08 Sep 2010 12:28:15 GMT

Hi
Does anybody have written a script or programme that creates a SAM or BAM file based on the following alignment information:

1) 454-reads
2) read qualtities
3) reference sequence(s)
4) read alignment information (e.g. read X aligns from (read-)position X_i to position X_j with reference Y, position Y_i to Y_j )

Many thanks for suggestions!

How to manage overlapping paired-end reads?

FredOnSeq — Wed, 08 Sep 2010 08:38:45 GMT

Hello,
This is my first time analysing RNA-Seq data; so I would appreciate your help with a couple of issues.
My data consist of 8 samples sequenced using Illumina�s standard paired-end RNA-seq protocol (which unfortunately at the time was NOT carried out in a strand-specific manner). The fragment size was 220bp which left around 100 bp after subtraction of adapters (2X58-60bp). Each pair (60 bp) in the dataset therefore has an approximately 20bp overlap.
I would like to use the bowtie/tophat/cufflinks pipeline and have a couple of questions regarding the analysis:
1. Although this post was very informative, I cannot decide which of these is a better strategy for analysing reads in my case:
a. Assembling the paired-end reads into 100bp single reads before�
b. Directly using the paired-end reads in tophat [This generates a negative (-20) inner distance between pairs but version 1.0.13 onwards seems to be able to handle this scenario (the �r option, I�m right???)].
2. Since the Illumina protocol was not strand-specific, is it a good idea to convert (correct word?) all the resulting mapping data or the sequence reads (after an initial round of mapping) so that it matches a single strand of the genome? I wonder if this strategy will help cufflinks better assemble and quantify the transcripts�
Thank you very much in advance for your help/suggestions/feedback�
Fred

variant annotation

Balat — Wed, 08 Sep 2010 04:12:12 GMT

Hi,
Most of the tools related to variant annotations seem to be related to humans. I work on a plant species where we only have the draft genome and not annotated yet. I have done an rna-seq experiment and I have found several snps using samtools. I got a gtf file by aligning reads against the genome sequence with bowtie and tophat. The gtf file only has transcirpt information and no orf or CDS information. Does any one have a script which takes the positions of the snps, annotations from a gtf file and the genome sequence in fasta format and predict if the snps are synonymous or non-synonymous?

SpliceMap 3.3.3.x released: Tool for spliced read alignment

john_mu — Tue, 07 Sep 2010 22:34:52 GMT

Hi Everyone,

SpliceMap 3.3.3.x has been release with many user requested features and stability/reliability enhancements.

Website: http://www.stanford.edu/group/wonglab/SpliceMap/

The following is a summary of the major changes.

High specificity in junction identification
Alignments of reads with variable length in both pairs (eg. after trimming)
Support for concatenated genome files
Bowtie index automatically built if not found
Read quality and names copied to SAM file
Option to run multiple chromosomes at same time
Option to change location of temp and output directories
and more... (see website)

A useful feature of SpliceMap I would just like to point out is that each junction is tagged with the number of non-redundant reads so that the reliability of each junction can be judged.

The following is some example output of junctions (not full alignment) from cisGenome browser.

Novel junctions are colored and possibly unreliable ones are colored more lightly.

Enjoy!

rna2map to blastparsed

ilyagoldin — Tue, 07 Sep 2010 22:12:19 GMT

Hello all. I've used rna2map to do an alignment, and now I want to feed the output to various tools that expect BLAST formatted input. How can I do that conversion? Thanks in advance.

Scripture score function

wuj — Tue, 07 Sep 2010 19:45:10 GMT

Hi, Everyone,

I am trying to use Scripture to do transcript abundance estimation, like Cufflinks does. I found there is a score function in the guide, so I simply use it with a transcription annotation database i have. I then use the RPKM column to calculate the expression of that transcript...But I after i compared the results with the outputs of cufflinks, they are not correlated....Did I used the wrong column, anyone has any suggestions? Thanks a lot!

http://www.broadinstitute.org/softwa...e/Score%20task

above is the link to the guide i mentioned.

UCSC for Arabidopsis

BioTalk — Tue, 07 Sep 2010 19:33:45 GMT

Hello All,

I would like to load the genome of Arabidopsis in UCSC browser using custom tracks. I am using the gff file available on the internet from ftp://ftp.arabidopsis.org/home/tair/...ransposons.gff

Does anyone know how to load genome for Arabidopsis in UCSC browser?

Thanks,

Bismark Bisulfite Aligner - Now supporting CpG, CHG and CHH context

fkrueger — Tue, 07 Sep 2010 15:33:03 GMT

I would like to announce that we have just released a new version of our Bisulfite-seq alignment and methylation calling tool Bismark. All associated files are available for free from http://www.bioinformatics.bbsrc.ac.uk/projects/.

As the most noticeable difference, Bismark does now further subdivide non-CpG context into CHG and CHH context, which will be especially interesting for researchers working on plant systems. The former characters 'C/c' in the methylation call has been replaced by:

CHG-context: X / x (methylated / unmethylated)
CHH-context: H / h (methylated / unmethylated)

In addition, I noticed that due to recent changes in the Bowtie source code, Bismark was producing lots of warnings 'best-first memory chunk exhaustion...') which was also mentioned here on SEQanswers. As suggested by Ben Langmead, the best way to counteract this problem is to increase the memory size for each bowtie thread, or mute bowtie. Thus, Bismark will now understand the additional option '--chunkmbs ' to change the memory from 64 (default)to any integer (I found that 256 or even 512 got rid of nearly all warnings). These errors were especially frequent in --best mode or for paired-end alignments. Bismark will now also understand the '--quiet' option to suppress memory chunk exhaustion (and other) warnings.

Some other minor fixes include:

- FastA files do no longer require the file ending ".fa".

- Fixed an issues so that Bismark will no longer tolerate chromosomes with
same name when reading the genome into memory.

- Fixed an issue with paired-end alignment reports.

- The methylation extractor will by default distinguish between cytosines in the three contexts CpG, CHG or CHH. If this is not needed, CHG and CHH context can be merged into 'non-CpG' context by specifying '--merge_non_CpG'.

- Due to the fundamental changes in v0.2.0 (CHG and CHH context methylation calls) the methylation extractor will now require that the Bismark mapping result file was generated with the same version of Bismark.

If you have any suggestion or comments I would like from you!

Best wishes,

Felix

Compile multi-individual SNP tables from bowtie/samtools mappings

KNS — Tue, 07 Sep 2010 14:07:38 GMT

I am sure someone must have solved this before, and I wonder if there is not an obviously easy way to do this:

I have GAIIx-RAD data of bowtie-reference-mapped individuals. The resulting individual SNP data (produced with samtools in the form of pileups) is easily merged in a single table, but one major problem remains:

The samtools pileup only reports on SNPs which are different wrt the reference sequence, but not on those SNPs which happen to be identical to the reference (i.e. those which show differences in other individuals).

Is there any easy way of discriminating the absence of data from DNA-sequence identity of any given SNP with the reference sequence?

Any pointers most welcome!:)

sequence capture visaulization

litali — Tue, 07 Sep 2010 12:14:59 GMT

Hello all
Does anybody know if there is a graphical solution to see the results of a sequence capture experiment. For example to see the targets and which regions were covered and to see the variations comparing to a reference file etc..
Thank you!!!

RNA expression with 454 data

litali — Tue, 07 Sep 2010 12:12:17 GMT

Is there any software to analyze 454 cDNA data to receive expression differences between samples?

hardware requirements to run blastX

litali — Tue, 07 Sep 2010 12:10:45 GMT

Hi all,
Can anybody give me an estimation what are the hardware requirements to have BlastX locally installed run in a reasonable pace and work ok with 454 data?
Many thanks!!!

GO annotation

Ramet — Tue, 07 Sep 2010 10:14:29 GMT

Hi all!

I have a problem with GO annotation of my BLAST results. So far I used Blast2GO, but for several reasons I am a bit unsatisfied with the programme.
In addition I need a good export function because I want to further analyse the data in MS Excel.

Can anyone name a sufficient programme which
- annotates GO terms for .xml files or swissprot accesion IDs (BLAST mustn't be done there!)
- possibility to only export level 2 GO terms
- possibility to export GO terms in tab-delimited form or equal for further analysis in Excel (information about category MF, BP, CC should be maintained!)

Thanks to everyone in advance!