Home / Browse / SAM tools / Mail / Archive

SAM tools

Email Archive: samtools-help (read-only)

2009:	_Jan	_Feb (26)	_Mar (37)	_Apr (20)	_May (20)	_Jun (81)	_Jul (65)	_Aug (77)	_Sep (85)	_Oct (86)	_Nov (233)	_Dec (131)
2010:	_Jan (201)	_Feb (274)	_Mar (257)	_Apr (208)	_May (99)	_Jun (223)	_Jul (131)	_Aug (139)	_Sep (195)	_Oct (133)	_Nov (196)	_Dec (184)
2011:	_Jan (156)	_Feb (181)	_Mar (246)	_Apr (131)	_May (172)	_Jun (157)	_Jul (121)	_Aug (200)	_Sep (199)	_Oct (114)	_Nov (189)	_Dec (68)
2012:	_Jan (130)	_Feb (180)	_Mar (232)	_Apr (173)	_May (149)	_Jun (105)	_Jul (107)	_Aug (69)	_Sep (59)	_Oct (49)	_Nov (43)	_Dec (56)
2013:	_Jan (69)	_Feb (58)	_Mar (79)	_Apr (90)	_May (61)	_Jun	_Jul	_Aug	_Sep	_Oct	_Nov	_Dec


S	M	T	W	T	F	S
	1	2	3	4	5	6
	(18)	(10)	(6)	(15)	(15)	(2)
7	8	9	10	11	12	13
	(12)	(25)	(21)	(7)	(12)
14	15	16	17	18	19	20
(3)		(15)	(4)	(7)	(15)	(2)
21	22	23	24	25	26	27
(2)	(20)	(9)	(10)	(21)	(9)	(7)
28
(7)

[Samtools-help] Subsetting BAM, at random.

From: Ivan Gregoretti <ivangreg@gm...> - 2010-03-03 01:34

Hi everybody,

How do you randomly subset a BAM file?

Lets say you have a BAM with 1000 records but you what 500 of them,
picked at random.

Can you set a random seed so that you can reproduce the subsetting?

Thank you,

Ivan


Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1592
Fax: 1-301-496-9878

Re: [Samtools-help] Subsetting BAM, at random.

From: Tom Blackwell <tblackw@um...> - 2010-03-03 01:48

man gawk | grep srand

On Tue, 2 Mar 2010, Ivan Gregoretti wrote:

> Hi everybody,
>
> How do you randomly subset a BAM file?
>
> Lets say you have a BAM with 1000 records but you what 500 of them,
> picked at random.
>
> Can you set a random seed so that you can reproduce the subsetting?
>
> Thank you,
>
> Ivan
>
>
> Ivan Gregoretti, PhD
> National Institute of Diabetes and Digestive and Kidney Diseases
> National Institutes of Health
> 5 Memorial Dr, Building 5, Room 205.
> Bethesda, MD 20892. USA.
> Phone: 1-301-496-1592
> Fax: 1-301-496-9878
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>
>
>

Re: [Samtools-help] Subsetting BAM, at random.

From: Shankar Ajay Subramanian <shankar.ajay@gm...> - 2010-03-03 02:01

Hi Ivan,

It seems like you want to sample a fraction (0.5) of the reads, in
which case, this might work:

samtools view <bamfile> | perl -e 'srand; while(<>){ print if
rand(1)<0.5 ; }' | <bamconversion>

Shankar

On Tue, Mar 2, 2010 at 8:34 PM, Ivan Gregoretti <ivangreg@...> wrote:
> Hi everybody,
>
> How do you randomly subset a BAM file?
>
> Lets say you have a BAM with 1000 records but you what 500 of them,
> picked at random.
>
> Can you set a random seed so that you can reproduce the subsetting?
>
> Thank you,
>
> Ivan
>
>
> Ivan Gregoretti, PhD
> National Institute of Diabetes and Digestive and Kidney Diseases
> National Institutes of Health
> 5 Memorial Dr, Building 5, Room 205.
> Bethesda, MD 20892. USA.
> Phone: 1-301-496-1592
> Fax: 1-301-496-9878
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>

Re: [Samtools-help] Subsetting BAM, at random.

From: Joseph Fass <joseph.fass@gm...> - 2010-03-03 02:09

Attachments: Message as HTML

or:

samtools view <bamfile> | perl -ne 'print if rand(1)<0.5' | <bamconversion>

(... which produces a different list every time ... though I'm not sure
about seed issues here, and exactly *how* random this is ...)

~Joe




On Tue, Mar 2, 2010 at 6:00 PM, Shankar Ajay Subramanian <
shankar.ajay@...> wrote:

> Hi Ivan,
>
> It seems like you want to sample a fraction (0.5) of the reads, in
> which case, this might work:
>
> samtools view <bamfile> | perl -e 'srand; while(<>){ print if
> rand(1)<0.5 ; }' | <bamconversion>
>
> Shankar
>
> On Tue, Mar 2, 2010 at 8:34 PM, Ivan Gregoretti <ivangreg@...>
> wrote:
> > Hi everybody,
> >
> > How do you randomly subset a BAM file?
> >
> > Lets say you have a BAM with 1000 records but you what 500 of them,
> > picked at random.
> >
> > Can you set a random seed so that you can reproduce the subsetting?
> >
> > Thank you,
> >
> > Ivan
> >
> >
> > Ivan Gregoretti, PhD
> > National Institute of Diabetes and Digestive and Kidney Diseases
> > National Institutes of Health
> > 5 Memorial Dr, Building 5, Room 205.
> > Bethesda, MD 20892. USA.
> > Phone: 1-301-496-1592
> > Fax: 1-301-496-9878
> >
> >
> ------------------------------------------------------------------------------
> > Download Intel&#174; Parallel Studio Eval
> > Try the new software tools for yourself. Speed compiling, find bugs
> > proactively, and fine-tune applications for parallel performance.
> > See why Intel Parallel Studio got high marks during beta.
> > http://p.sf.net/sfu/intel-sw-dev
> > _______________________________________________
> > Samtools-help mailing list
> > Samtools-help@...
> > https://lists.sourceforge.net/lists/listinfo/samtools-help
> >
>
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>



-- 
Joseph Fass
Bioinformatics Programmer
UC Davis Bioinformatics Core
joseph.fass -at- gmail.com (professional)
970.227.5928 (c)  ||  530.752.2698 (w)

Re: [Samtools-help] Subsetting BAM, at random.

From: Ivan Gregoretti <ivangreg@gm...> - 2010-03-03 02:36

Hi Shankar and Joe,

Actually, you have illustrated my point very well:

piping through perl does not solve the problem of the BAM header. I
want to recover the header of the original BAM without calling the
ref.fa.

In Shankar's example, this would be

samtools view fullset.bam | perl -e 'srand; while(<>){ print if
rand(1)<0.5 ; }' | samtools view -t ref.fa -bS - > subset.bam

I was wondering if there was any samtools subsetting functionality to
go from BAM to BAM without piping through an external program. Well,
apparently, there is not. Too bad; I would expect it a very common
operation.

Thank you for your examples. Keep saming and baming.

Ivan
PS: also, keep perling!


Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1592
Fax: 1-301-496-9878



On Tue, Mar 2, 2010 at 9:08 PM, Joseph Fass <joseph.fass@...> wrote:
> or:
>
> samtools view <bamfile> | perl -ne 'print if rand(1)<0.5' | <bamconversion>
>
> (... which produces a different list every time ... though I'm not sure
> about seed issues here, and exactly how random this is ...)
>
> ~Joe
>
>
>
>
> On Tue, Mar 2, 2010 at 6:00 PM, Shankar Ajay Subramanian
> <shankar.ajay@...> wrote:
>>
>> Hi Ivan,
>>
>> It seems like you want to sample a fraction (0.5) of the reads, in
>> which case, this might work:
>>
>> samtools view <bamfile> | perl -e 'srand; while(<>){ print if
>> rand(1)<0.5 ; }' | <bamconversion>
>>
>> Shankar
>>
>> On Tue, Mar 2, 2010 at 8:34 PM, Ivan Gregoretti <ivangreg@...>
>> wrote:
>> > Hi everybody,
>> >
>> > How do you randomly subset a BAM file?
>> >
>> > Lets say you have a BAM with 1000 records but you what 500 of them,
>> > picked at random.
>> >
>> > Can you set a random seed so that you can reproduce the subsetting?
>> >
>> > Thank you,
>> >
>> > Ivan
>> >
>> >
>> > Ivan Gregoretti, PhD
>> > National Institute of Diabetes and Digestive and Kidney Diseases
>> > National Institutes of Health
>> > 5 Memorial Dr, Building 5, Room 205.
>> > Bethesda, MD 20892. USA.
>> > Phone: 1-301-496-1592
>> > Fax: 1-301-496-9878
>> >
>> >
>> > ------------------------------------------------------------------------------
>> > Download Intel&#174; Parallel Studio Eval
>> > Try the new software tools for yourself. Speed compiling, find bugs
>> > proactively, and fine-tune applications for parallel performance.
>> > See why Intel Parallel Studio got high marks during beta.
>> > http://p.sf.net/sfu/intel-sw-dev
>> > _______________________________________________
>> > Samtools-help mailing list
>> > Samtools-help@...
>> > https://lists.sourceforge.net/lists/listinfo/samtools-help
>> >
>>
>>
>> ------------------------------------------------------------------------------
>> Download Intel&#174; Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Samtools-help mailing list
>> Samtools-help@...
>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>
>
>
> --
> Joseph Fass
> Bioinformatics Programmer
> UC Davis Bioinformatics Core
> joseph.fass -at- gmail.com (professional)
> 970.227.5928 (c)  ||  530.752.2698 (w)
>

Re: [Samtools-help] Subsetting BAM, at random.

From: Alec Wysoker <alecw@br...> - 2010-03-03 03:01

Hi Ivan,

You can either write a Java program using the Picard SAM library, or you 
can use samtools, convert the bam to sam, pass through the header, 
randomly select lines using some kind of text-processing method, and 
convert back into bam.

-Alec

Ivan Gregoretti wrote:
> Hi everybody,
>
> How do you randomly subset a BAM file?
>
> Lets say you have a BAM with 1000 records but you what 500 of them,
> picked at random.
>
> Can you set a random seed so that you can reproduce the subsetting?
>
> Thank you,
>
> Ivan
>
>
> Ivan Gregoretti, PhD
> National Institute of Diabetes and Digestive and Kidney Diseases
> National Institutes of Health
> 5 Memorial Dr, Building 5, Room 205.
> Bethesda, MD 20892. USA.
> Phone: 1-301-496-1592
> Fax: 1-301-496-9878
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>

Re: [Samtools-help] Subsetting BAM, at random.

From: Mark A. DePristo <depristo@br...> - 2010-03-03 03:33

You can do this with the GATK using PrintReads and the engine's  
downsampling arguments:

	http://www.broadinstitute.org/gsa/wiki/index.php/The_Genome_Analysis_Toolkit

Such as:

	java -jar GenomeAnalysisTK.jar -I in.bam -T PrintReads -- 
downsample_to_fraction 0.5 -o out.bam

Best,

Mark

On Mar 2, 2010, at 9:35 PM, Ivan Gregoretti wrote:

> Hi Shankar and Joe,
>
> Actually, you have illustrated my point very well:
>
> piping through perl does not solve the problem of the BAM header. I
> want to recover the header of the original BAM without calling the
> ref.fa.
>
> In Shankar's example, this would be
>
> samtools view fullset.bam | perl -e 'srand; while(<>){ print if
> rand(1)<0.5 ; }' | samtools view -t ref.fa -bS - > subset.bam
>
> I was wondering if there was any samtools subsetting functionality to
> go from BAM to BAM without piping through an external program. Well,
> apparently, there is not. Too bad; I would expect it a very common
> operation.
>
> Thank you for your examples. Keep saming and baming.
>
> Ivan
> PS: also, keep perling!
>
>
> Ivan Gregoretti, PhD
> National Institute of Diabetes and Digestive and Kidney Diseases
> National Institutes of Health
> 5 Memorial Dr, Building 5, Room 205.
> Bethesda, MD 20892. USA.
> Phone: 1-301-496-1592
> Fax: 1-301-496-9878
>
>
>
> On Tue, Mar 2, 2010 at 9:08 PM, Joseph Fass <joseph.fass@...>  
> wrote:
>> or:
>>
>> samtools view <bamfile> | perl -ne 'print if rand(1)<0.5' |  
>> <bamconversion>
>>
>> (... which produces a different list every time ... though I'm not  
>> sure
>> about seed issues here, and exactly how random this is ...)
>>
>> ~Joe
>>
>>
>>
>>
>> On Tue, Mar 2, 2010 at 6:00 PM, Shankar Ajay Subramanian
>> <shankar.ajay@...> wrote:
>>>
>>> Hi Ivan,
>>>
>>> It seems like you want to sample a fraction (0.5) of the reads, in
>>> which case, this might work:
>>>
>>> samtools view <bamfile> | perl -e 'srand; while(<>){ print if
>>> rand(1)<0.5 ; }' | <bamconversion>
>>>
>>> Shankar
>>>
>>> On Tue, Mar 2, 2010 at 8:34 PM, Ivan Gregoretti <ivangreg@...>
>>> wrote:
>>>> Hi everybody,
>>>>
>>>> How do you randomly subset a BAM file?
>>>>
>>>> Lets say you have a BAM with 1000 records but you what 500 of them,
>>>> picked at random.
>>>>
>>>> Can you set a random seed so that you can reproduce the subsetting?
>>>>
>>>> Thank you,
>>>>
>>>> Ivan
>>>>
>>>>
>>>> Ivan Gregoretti, PhD
>>>> National Institute of Diabetes and Digestive and Kidney Diseases
>>>> National Institutes of Health
>>>> 5 Memorial Dr, Building 5, Room 205.
>>>> Bethesda, MD 20892. USA.
>>>> Phone: 1-301-496-1592
>>>> Fax: 1-301-496-9878
>>>>
>>>>
>>>> ------------------------------------------------------------------------------
>>>> Download Intel&#174; Parallel Studio Eval
>>>> Try the new software tools for yourself. Speed compiling, find bugs
>>>> proactively, and fine-tune applications for parallel performance.
>>>> See why Intel Parallel Studio got high marks during beta.
>>>> http://p.sf.net/sfu/intel-sw-dev
>>>> _______________________________________________
>>>> Samtools-help mailing list
>>>> Samtools-help@...
>>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>>>
>>>
>>>
>>> ------------------------------------------------------------------------------
>>> Download Intel&#174; Parallel Studio Eval
>>> Try the new software tools for yourself. Speed compiling, find bugs
>>> proactively, and fine-tune applications for parallel performance.
>>> See why Intel Parallel Studio got high marks during beta.
>>> http://p.sf.net/sfu/intel-sw-dev
>>> _______________________________________________
>>> Samtools-help mailing list
>>> Samtools-help@...
>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>
>>
>>
>> --
>> Joseph Fass
>> Bioinformatics Programmer
>> UC Davis Bioinformatics Core
>> joseph.fass -at- gmail.com (professional)
>> 970.227.5928 (c)  ||  530.752.2698 (w)
>>
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help

Mark A. DePristo, Ph.D.
Manager, Medical and Population Genetics Analysis
Broad Institute of Harvard and MIT
depristo@...
mark@...

Re: [Samtools-help] Subsetting BAM, at random.

From: Ivan Gregoretti <ivangreg@gm...> - 2010-03-03 07:50

Hi Mark,

That looks great. Does it work on non-indexed bam files?

Thanks,

Ivan

Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1592
Fax: 1-301-496-9878



On Tue, Mar 2, 2010 at 10:33 PM, Mark A. DePristo
<depristo@...> wrote:
> You can do this with the GATK using PrintReads and the engine's downsampling
> arguments:
>
>
>  http://www.broadinstitute.org/gsa/wiki/index.php/The_Genome_Analysis_Toolkit
>
> Such as:
>
>        java -jar GenomeAnalysisTK.jar -I in.bam -T PrintReads
> --downsample_to_fraction 0.5 -o out.bam
>
> Best,
>
> Mark
>
> On Mar 2, 2010, at 9:35 PM, Ivan Gregoretti wrote:
>
>> Hi Shankar and Joe,
>>
>> Actually, you have illustrated my point very well:
>>
>> piping through perl does not solve the problem of the BAM header. I
>> want to recover the header of the original BAM without calling the
>> ref.fa.
>>
>> In Shankar's example, this would be
>>
>> samtools view fullset.bam | perl -e 'srand; while(<>){ print if
>> rand(1)<0.5 ; }' | samtools view -t ref.fa -bS - > subset.bam
>>
>> I was wondering if there was any samtools subsetting functionality to
>> go from BAM to BAM without piping through an external program. Well,
>> apparently, there is not. Too bad; I would expect it a very common
>> operation.
>>
>> Thank you for your examples. Keep saming and baming.
>>
>> Ivan
>> PS: also, keep perling!
>>
>>
>> Ivan Gregoretti, PhD
>> National Institute of Diabetes and Digestive and Kidney Diseases
>> National Institutes of Health
>> 5 Memorial Dr, Building 5, Room 205.
>> Bethesda, MD 20892. USA.
>> Phone: 1-301-496-1592
>> Fax: 1-301-496-9878
>>
>>
>>
>> On Tue, Mar 2, 2010 at 9:08 PM, Joseph Fass <joseph.fass@...> wrote:
>>>
>>> or:
>>>
>>> samtools view <bamfile> | perl -ne 'print if rand(1)<0.5' |
>>> <bamconversion>
>>>
>>> (... which produces a different list every time ... though I'm not sure
>>> about seed issues here, and exactly how random this is ...)
>>>
>>> ~Joe
>>>
>>>
>>>
>>>
>>> On Tue, Mar 2, 2010 at 6:00 PM, Shankar Ajay Subramanian
>>> <shankar.ajay@...> wrote:
>>>>
>>>> Hi Ivan,
>>>>
>>>> It seems like you want to sample a fraction (0.5) of the reads, in
>>>> which case, this might work:
>>>>
>>>> samtools view <bamfile> | perl -e 'srand; while(<>){ print if
>>>> rand(1)<0.5 ; }' | <bamconversion>
>>>>
>>>> Shankar
>>>>
>>>> On Tue, Mar 2, 2010 at 8:34 PM, Ivan Gregoretti <ivangreg@...>
>>>> wrote:
>>>>>
>>>>> Hi everybody,
>>>>>
>>>>> How do you randomly subset a BAM file?
>>>>>
>>>>> Lets say you have a BAM with 1000 records but you what 500 of them,
>>>>> picked at random.
>>>>>
>>>>> Can you set a random seed so that you can reproduce the subsetting?
>>>>>
>>>>> Thank you,
>>>>>
>>>>> Ivan
>>>>>
>>>>>
>>>>> Ivan Gregoretti, PhD
>>>>> National Institute of Diabetes and Digestive and Kidney Diseases
>>>>> National Institutes of Health
>>>>> 5 Memorial Dr, Building 5, Room 205.
>>>>> Bethesda, MD 20892. USA.
>>>>> Phone: 1-301-496-1592
>>>>> Fax: 1-301-496-9878
>>>>>
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>> Download Intel&#174; Parallel Studio Eval
>>>>> Try the new software tools for yourself. Speed compiling, find bugs
>>>>> proactively, and fine-tune applications for parallel performance.
>>>>> See why Intel Parallel Studio got high marks during beta.
>>>>> http://p.sf.net/sfu/intel-sw-dev
>>>>> _______________________________________________
>>>>> Samtools-help mailing list
>>>>> Samtools-help@...
>>>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>>>>
>>>>
>>>>
>>>>
>>>> ------------------------------------------------------------------------------
>>>> Download Intel&#174; Parallel Studio Eval
>>>> Try the new software tools for yourself. Speed compiling, find bugs
>>>> proactively, and fine-tune applications for parallel performance.
>>>> See why Intel Parallel Studio got high marks during beta.
>>>> http://p.sf.net/sfu/intel-sw-dev
>>>> _______________________________________________
>>>> Samtools-help mailing list
>>>> Samtools-help@...
>>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>>
>>>
>>>
>>> --
>>> Joseph Fass
>>> Bioinformatics Programmer
>>> UC Davis Bioinformatics Core
>>> joseph.fass -at- gmail.com (professional)
>>> 970.227.5928 (c)  ||  530.752.2698 (w)
>>>
>>
>>
>> ------------------------------------------------------------------------------
>> Download Intel&#174; Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Samtools-help mailing list
>> Samtools-help@...
>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>
> Mark A. DePristo, Ph.D.
> Manager, Medical and Population Genetics Analysis
> Broad Institute of Harvard and MIT
> depristo@...
> mark@...
>
>
>
>
>
>
>
>
>
>
>
>
>
>

[Samtools-help] tview displaying N's for reference when no mapped reads in view

From: Joseph Fass <joseph.fass@gm...> - 2010-03-02 19:24

Attachments: Message as HTML

I'm seeing something odd with tview, with bam files produced from SOLiD
reads with bwa-0.5.7 and tview'd with samtools 0.1.7 ... it seems that when
no mapped reads are in view, all reference bases change to N's ... so:

         451       461       471       481       491       501
511       521       531
ttttcttcttttctttttcttatttttctttttacagttttagttttcttttatttttgaatagaaattagctttataaaatatagtaacaccacctaaattggtt

.........

,,,,,,,,,

... becomes the following after shifting a to the left so that the read is
off-screen:

        441       451       461       471       481       491
501       511       521
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN


... the N's change back to ATCG's if only one mapped base is visible.

The situation from the left is different ... the reference bases to the
right of the last mapped base are always N's (perhaps only when the unmapped
range is longer than one screen? -- I'm not sure).

~Joe





-- 
Joseph Fass
Bioinformatics Programmer
UC Davis Bioinformatics Core
joseph.fass -at- gmail.com (professional)
970.227.5928 (c)  ||  530.752.2698 (w)

Re: [Samtools-help] tview displaying N's for reference when no mapped reads in view

From: Joseph Fass <joseph.fass@gm...> - 2010-03-03 02:23

Attachments: Message as HTML

I've now seen this behavior with 454 reads ... so it seems to be independent
of bwa ...

~Joe




On Tue, Mar 2, 2010 at 11:24 AM, Joseph Fass <joseph.fass@...> wrote:

> I'm seeing something odd with tview, with bam files produced from SOLiD
> reads with bwa-0.5.7 and tview'd with samtools 0.1.7 ... it seems that when
> no mapped reads are in view, all reference bases change to N's ... so:
>
>          451       461       471       481       491       501
> 511       521       531
> ttttcttcttttctttttcttatttttctttttacagttttagttttcttttatttttgaatagaaattagctttataaaatatagtaacaccacctaaattggtt
>
> .........
>
> ,,,,,,,,,
>
> ... becomes the following after shifting a to the left so that the read is
> off-screen:
>
>         441       451       461       471       481       491
> 501       511       521
> NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>
>
> ... the N's change back to ATCG's if only one mapped base is visible.
>
> The situation from the left is different ... the reference bases to the
> right of the last mapped base are always N's (perhaps only when the unmapped
> range is longer than one screen? -- I'm not sure).
>
> ~Joe
>
>
>
>
>
> --
> Joseph Fass
> Bioinformatics Programmer
> UC Davis Bioinformatics Core
> joseph.fass -at- gmail.com (professional)
> 970.227.5928 (c)  ||  530.752.2698 (w)
>



-- 
Joseph Fass
Bioinformatics Programmer
UC Davis Bioinformatics Core
joseph.fass -at- gmail.com (professional)
970.227.5928 (c)  ||  530.752.2698 (w)

Re: [Samtools-help] tview displaying N's for reference when no mapped reads in view

From: Kevin Lam <aboulia@gm...> - 2010-03-03 03:17

Attachments: Message as HTML

I have the same problem as well.
but with
SOLiD v3.5 transcriptome reads
bwa-0.5.6
samtools-0.1.7_x86_64-linux

when I use wgsim from samtools to generate the data though all is fine. Have
you tried that?

I am going to try bowtie to see if it's a BWA problem or SOLiD problem.
:(

Kevin

On Wed, Mar 3, 2010 at 3:24 AM, Joseph Fass <joseph.fass@...> wrote:

> I'm seeing something odd with tview, with bam files produced from SOLiD
> reads with bwa-0.5.7 and tview'd with samtools 0.1.7 ... it seems that when
> no mapped reads are in view, all reference bases change to N's ... so:
>
>          451       461       471       481       491       501
> 511       521       531
> ttttcttcttttctttttcttatttttctttttacagttttagttttcttttatttttgaatagaaattagctttataaaatatagtaacaccacctaaattggtt
>
> .........
>
> ,,,,,,,,,
>
> ... becomes the following after shifting a to the left so that the read is
> off-screen:
>
>         441       451       461       471       481       491
> 501       511       521
> NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>
>
> ... the N's change back to ATCG's if only one mapped base is visible.
>
> The situation from the left is different ... the reference bases to the
> right of the last mapped base are always N's (perhaps only when the unmapped
> range is longer than one screen? -- I'm not sure).
>
> ~Joe
>
>
>
>
>
> --
> Joseph Fass
> Bioinformatics Programmer
> UC Davis Bioinformatics Core
> joseph.fass -at- gmail.com (professional)
> 970.227.5928 (c)  ||  530.752.2698 (w)
>
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>
>

[Samtools-help] Bio::DB::Sam question - filtering coverage data

From: Keiran Raine <kr2@sa...> - 2010-03-02 14:50

Attachments: Message as HTML

Hi,

Does anyone know if it is possible to filter coverage data based on  
the flag field?

What I want is the coverage over a region where the 'mate' is mapped.   
The code below appears to work, but I want to be sure that it is not  
just ignoring the '-filter' component as it's not entirely clear in  
the documentation.


	my $base_count = 0;
	my $depth = 0;
	my ($coverage) = $sam->features(-type=>'coverage',-seq_id=>$chr, - 
start=>$from, -end=>$to,
									-filter=> sub {my $flag = shift->flag; return ($flag | 8) !=  
$flag});;
	if($coverage) {
		foreach my $curr_depth(@{$coverage->coverage}) {
			$depth+=$curr_depth;
			$base_count++;
		}
	}
	if($base_count > 0) {
		$depth = int ($depth/$base_count);
	}

Keiran Raine
Senior Computer Biologist
The Cancer Genome Project
Ext: 2100
kr2@...

The Wellcome Trust Sanger Institute is operated by Genome Research
Limited, a charity registered in England with number 1021457 and a
company registered in England with number 2742969, whose registered
office is 215 Euston Road, London, NW1 2BE.






-- 
 The Wellcome Trust Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE.

[Samtools-help] Picard release 1.15

From: Alec Wysoker <alecw@br...> - 2010-03-01 17:56

Release notes for Picard 1.15
1 March 2010

- Add code to write LiftOver Chain file.
- Add validation to LiftOver.
- In CommandLineParser, allow a special symbol to be used in parameter 
names to force the option to be used even if coming from the params file 
and would usually be overridden by the command line.
- Remove final from various public attribute-related SAMRecord methods.  
I'm not sure why these were final to begin with.
- Added support for clipping to SamToFastq.
- Refactored cigar-clipping into a separate iterator for use outside of 
command-line program CleanSam.
- Add utility to reconstruct reference bases from read with cigar and MD 
string.  Not tested exhaustively, so please report any problems.
- Warn if paired read is not marked as first or second of pair.
- Add SamLocusIterator -- class for iterating over a SAM file by locus, 
and getting pile of reads and offsets that align with locus.
- Add getter and setter for predicted median insert size to 
SAMReadGroupRecord.
- Modified ValidateSamFile to return a non-zero code when validation fails.
- Change ValidateSamFile so that invalid value for @RG PI tag triggers 
error message instead of exception.

-Alec

Re: [Samtools-help] Picard release 1.15

From: Ivan Gregoretti <ivangreg@gm...> - 2010-03-01 19:38

Hi again,

This time I bring up an issue of bigger consequence that I just
spotted in this release. It has to do with boundary conditions on BAM
files. This is the example:

I have a SAM record that goes like this

zcat s_5.sam.gz | perl -wlne'if($_=~/HWI-EAS179_1:5:71:1726:6852*/){print}'
HWI-EAS179_1:5:71:1726:6852	0	chrM	16266	118	36M	*	0	0	TAATCATACTCTATTACGCAATAAACATTAACAAGT	?>:;77/8<7A=BAABDDA=BB@...	XD:Z:34NN	SM:i:118

As you can see, it is 36 bp record that starts at position 16266 and
it matches the forward strand of chrM.

Now, the corresponding BAM file was successfully created like this

samtools view -t mouse.mm9.genome -bS s_5.sam.gz > s_5.bam

where an excerpt of 'cat mouse.mm9.genome' goes like this
...
chrY	15902555
chrM	16299
chr13_random	400311
...

Notice that 36 + 16266 = 16302. Ugh.
ViewSAM.jar complains that this records is out of range and terminates parsing.

Now, mitochondrial DNA is linear!

Ivan


Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1592
Fax: 1-301-496-9878



On Mon, Mar 1, 2010 at 12:55 PM, Alec Wysoker <alecw@...> wrote:
> Release notes for Picard 1.15
> 1 March 2010
>
> - Add code to write LiftOver Chain file.
> - Add validation to LiftOver.
> - In CommandLineParser, allow a special symbol to be used in parameter
> names to force the option to be used even if coming from the params file
> and would usually be overridden by the command line.
> - Remove final from various public attribute-related SAMRecord methods.
> I'm not sure why these were final to begin with.
> - Added support for clipping to SamToFastq.
> - Refactored cigar-clipping into a separate iterator for use outside of
> command-line program CleanSam.
> - Add utility to reconstruct reference bases from read with cigar and MD
> string.  Not tested exhaustively, so please report any problems.
> - Warn if paired read is not marked as first or second of pair.
> - Add SamLocusIterator -- class for iterating over a SAM file by locus,
> and getting pile of reads and offsets that align with locus.
> - Add getter and setter for predicted median insert size to
> SAMReadGroupRecord.
> - Modified ValidateSamFile to return a non-zero code when validation fails.
> - Change ValidateSamFile so that invalid value for @RG PI tag triggers
> error message instead of exception.
>
> -Alec
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>

Re: [Samtools-help] [Samtools-announce] Picard release 1.15

From: Alec Wysoker <alecw@br...> - 2010-03-01 22:09

Hi Klaudia,

 From the help text for SamToFastq:

CLIPPING_ATTRIBUTE=String
CLIP_ATTR=String              The attribute that stores the position at 
which the SAM record should be clipped. Default
                              value: null.

CLIPPING_ACTION=String
CLIP_ACT=String               The action that should be taken with 
clipped reads: 'X' means the reads and qualities
                              should be trimmed at the clipped position; 
'N' means the bases should be changed to Ns in
                              the clipped region; and any integer means 
that the base qualities should be set to that
                              value in the clipped region. Default 
value: null.


-Alec

Klaudia Walter wrote:
> Hi Alec,
>
> What do you mean by
>
> - Added support for clipping to SamToFastq.
>
> Thanks,
> Klaudia
>
>
> On 1 Mar 2010, at 17:55, Alec Wysoker wrote:
>
>> Release notes for Picard 1.15
>> 1 March 2010
>>
>> - Add code to write LiftOver Chain file.
>> - Add validation to LiftOver.
>> - In CommandLineParser, allow a special symbol to be used in parameter
>> names to force the option to be used even if coming from the params file
>> and would usually be overridden by the command line.
>> - Remove final from various public attribute-related SAMRecord methods.
>> I'm not sure why these were final to begin with.
>> - Added support for clipping to SamToFastq.
>> - Refactored cigar-clipping into a separate iterator for use outside of
>> command-line program CleanSam.
>> - Add utility to reconstruct reference bases from read with cigar and MD
>> string.  Not tested exhaustively, so please report any problems.
>> - Warn if paired read is not marked as first or second of pair.
>> - Add SamLocusIterator -- class for iterating over a SAM file by locus,
>> and getting pile of reads and offsets that align with locus.
>> - Add getter and setter for predicted median insert size to
>> SAMReadGroupRecord.
>> - Modified ValidateSamFile to return a non-zero code when validation 
>> fails.
>> - Change ValidateSamFile so that invalid value for @RG PI tag triggers
>> error message instead of exception.
>>
>> -Alec
>>
>> ------------------------------------------------------------------------------ 
>>
>> Download Intel&#174; Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Samtools-announce mailing list
>> Samtools-announce@...
>> https://lists.sourceforge.net/lists/listinfo/samtools-announce
>
>
>

Re: [Samtools-help] Picard release 1.15

From: Alec Wysoker <alecw@br...> - 2010-03-01 22:14

Hi Ivan,

Try passing VALIDATION_STRINGENCY=LENIENT or 
VALIDATION_STRINGENCY=SILENT to ViewSam.

-Alec

Ivan Gregoretti wrote:
> Hi again,
>
> This time I bring up an issue of bigger consequence that I just
> spotted in this release. It has to do with boundary conditions on BAM
> files. This is the example:
>
> I have a SAM record that goes like this
>
> zcat s_5.sam.gz | perl -wlne'if($_=~/HWI-EAS179_1:5:71:1726:6852*/){print}'
> HWI-EAS179_1:5:71:1726:6852	0	chrM	16266	118	36M	*	0	0	TAATCATACTCTATTACGCAATAAACATTAACAAGT	?>:;77/8<7A=BAABDDA=BB@...	XD:Z:34NN	SM:i:118
>
> As you can see, it is 36 bp record that starts at position 16266 and
> it matches the forward strand of chrM.
>
> Now, the corresponding BAM file was successfully created like this
>
> samtools view -t mouse.mm9.genome -bS s_5.sam.gz > s_5.bam
>
> where an excerpt of 'cat mouse.mm9.genome' goes like this
> ...
> chrY	15902555
> chrM	16299
> chr13_random	400311
> ...
>
> Notice that 36 + 16266 = 16302. Ugh.
> ViewSAM.jar complains that this records is out of range and terminates parsing.
>
> Now, mitochondrial DNA is linear!
>
> Ivan
>
>
> Ivan Gregoretti, PhD
> National Institute of Diabetes and Digestive and Kidney Diseases
> National Institutes of Health
> 5 Memorial Dr, Building 5, Room 205.
> Bethesda, MD 20892. USA.
> Phone: 1-301-496-1592
> Fax: 1-301-496-9878
>
>
>
> On Mon, Mar 1, 2010 at 12:55 PM, Alec Wysoker <alecw@...> wrote:
>   
>> Release notes for Picard 1.15
>> 1 March 2010
>>
>> - Add code to write LiftOver Chain file.
>> - Add validation to LiftOver.
>> - In CommandLineParser, allow a special symbol to be used in parameter
>> names to force the option to be used even if coming from the params file
>> and would usually be overridden by the command line.
>> - Remove final from various public attribute-related SAMRecord methods.
>> I'm not sure why these were final to begin with.
>> - Added support for clipping to SamToFastq.
>> - Refactored cigar-clipping into a separate iterator for use outside of
>> command-line program CleanSam.
>> - Add utility to reconstruct reference bases from read with cigar and MD
>> string.  Not tested exhaustively, so please report any problems.
>> - Warn if paired read is not marked as first or second of pair.
>> - Add SamLocusIterator -- class for iterating over a SAM file by locus,
>> and getting pile of reads and offsets that align with locus.
>> - Add getter and setter for predicted median insert size to
>> SAMReadGroupRecord.
>> - Modified ValidateSamFile to return a non-zero code when validation fails.
>> - Change ValidateSamFile so that invalid value for @RG PI tag triggers
>> error message instead of exception.
>>
>> -Alec
>>
>> ------------------------------------------------------------------------------
>> Download Intel&#174; Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Samtools-help mailing list
>> Samtools-help@...
>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>
>>

Re: [Samtools-help] Picard release 1.15

From: Ivan Gregoretti <ivangreg@gm...> - 2010-03-02 02:45

Thanks Alec. I'll try that.

By the way, when in my post I said that
"mitochondrial DNA is linear",
I actually meant that it is circular. I see that you got it anyway.

Cheers,

Ivan

Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1592
Fax: 1-301-496-9878



On Mon, Mar 1, 2010 at 5:14 PM, Alec Wysoker <alecw@...> wrote:
> Hi Ivan,
>
> Try passing VALIDATION_STRINGENCY=LENIENT or VALIDATION_STRINGENCY=SILENT to
> ViewSam.
>
> -Alec
>
> Ivan Gregoretti wrote:
>>
>> Hi again,
>>
>> This time I bring up an issue of bigger consequence that I just
>> spotted in this release. It has to do with boundary conditions on BAM
>> files. This is the example:
>>
>> I have a SAM record that goes like this
>>
>> zcat s_5.sam.gz | perl
>> -wlne'if($_=~/HWI-EAS179_1:5:71:1726:6852*/){print}'
>> HWI-EAS179_1:5:71:1726:6852     0       chrM    16266   118     36M     *
>>       0       0       TAATCATACTCTATTACGCAATAAACATTAACAAGT
>>  ?>:;77/8<7A=BAABDDA=BB@...    XD:Z:34NN       SM:i:118
>>
>> As you can see, it is 36 bp record that starts at position 16266 and
>> it matches the forward strand of chrM.
>>
>> Now, the corresponding BAM file was successfully created like this
>>
>> samtools view -t mouse.mm9.genome -bS s_5.sam.gz > s_5.bam
>>
>> where an excerpt of 'cat mouse.mm9.genome' goes like this
>> ...
>> chrY    15902555
>> chrM    16299
>> chr13_random    400311
>> ...
>>
>> Notice that 36 + 16266 = 16302. Ugh.
>> ViewSAM.jar complains that this records is out of range and terminates
>> parsing.
>>
>> Now, mitochondrial DNA is linear!
>>
>> Ivan
>>
>>
>> Ivan Gregoretti, PhD
>> National Institute of Diabetes and Digestive and Kidney Diseases
>> National Institutes of Health
>> 5 Memorial Dr, Building 5, Room 205.
>> Bethesda, MD 20892. USA.
>> Phone: 1-301-496-1592
>> Fax: 1-301-496-9878
>>
>>
>>
>> On Mon, Mar 1, 2010 at 12:55 PM, Alec Wysoker <alecw@...>
>> wrote:
>>
>>>
>>> Release notes for Picard 1.15
>>> 1 March 2010
>>>
>>> - Add code to write LiftOver Chain file.
>>> - Add validation to LiftOver.
>>> - In CommandLineParser, allow a special symbol to be used in parameter
>>> names to force the option to be used even if coming from the params file
>>> and would usually be overridden by the command line.
>>> - Remove final from various public attribute-related SAMRecord methods.
>>> I'm not sure why these were final to begin with.
>>> - Added support for clipping to SamToFastq.
>>> - Refactored cigar-clipping into a separate iterator for use outside of
>>> command-line program CleanSam.
>>> - Add utility to reconstruct reference bases from read with cigar and MD
>>> string.  Not tested exhaustively, so please report any problems.
>>> - Warn if paired read is not marked as first or second of pair.
>>> - Add SamLocusIterator -- class for iterating over a SAM file by locus,
>>> and getting pile of reads and offsets that align with locus.
>>> - Add getter and setter for predicted median insert size to
>>> SAMReadGroupRecord.
>>> - Modified ValidateSamFile to return a non-zero code when validation
>>> fails.
>>> - Change ValidateSamFile so that invalid value for @RG PI tag triggers
>>> error message instead of exception.
>>>
>>> -Alec
>>>
>>>
>>> ------------------------------------------------------------------------------
>>> Download Intel&#174; Parallel Studio Eval
>>> Try the new software tools for yourself. Speed compiling, find bugs
>>> proactively, and fine-tune applications for parallel performance.
>>> See why Intel Parallel Studio got high marks during beta.
>>> http://p.sf.net/sfu/intel-sw-dev
>>> _______________________________________________
>>> Samtools-help mailing list
>>> Samtools-help@...
>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>>
>>>
>

Re: [Samtools-help] Picard release 1.15

From: Keiran Raine <kr2@sa...> - 2010-03-02 08:30

Hi,

What you are seeing is actually a bug/feature from BWA (which should  
not occur in the latest release).

Your read is actually mapping to the last 33 bases of chrM and the  
first 3 bases of chr13_random.  Under the latest release this read  
would be marked as unmapped.

Regards,

Keiran Raine
Senior Computer Biologist
The Cancer Genome Project
Ext: 2100
kr2@...

The Wellcome Trust Sanger Institute is operated by Genome Research
Limited, a charity registered in England with number 1021457 and a
company registered in England with number 2742969, whose registered
office is 215 Euston Road, London, NW1 2BE.



On 1 Mar 2010, at 19:38, Ivan Gregoretti wrote:

> Hi again,
>
> This time I bring up an issue of bigger consequence that I just
> spotted in this release. It has to do with boundary conditions on BAM
> files. This is the example:
>
> I have a SAM record that goes like this
>
> zcat s_5.sam.gz | perl -wlne'if($_=~/HWI-EAS179_1:5:71:1726:6852*/) 
> {print}'
> HWI-EAS179_1:5:71:1726:6852	0	chrM	16266	118	36M	*	0	0	 
> TAATCATACTCTATTACGCAATAAACATTAACAAGT	?>:;77/8<7A=BAABDDA=BB@...? 
> B<DB	XD:Z:34NN	SM:i:118
>
> As you can see, it is 36 bp record that starts at position 16266 and
> it matches the forward strand of chrM.
>
> Now, the corresponding BAM file was successfully created like this
>
> samtools view -t mouse.mm9.genome -bS s_5.sam.gz > s_5.bam
>
> where an excerpt of 'cat mouse.mm9.genome' goes like this
> ...
> chrY	15902555
> chrM	16299
> chr13_random	400311
> ...
>
> Notice that 36 + 16266 = 16302. Ugh.
> ViewSAM.jar complains that this records is out of range and  
> terminates parsing.
>
> Now, mitochondrial DNA is linear!
>
> Ivan
>
>
> Ivan Gregoretti, PhD
> National Institute of Diabetes and Digestive and Kidney Diseases
> National Institutes of Health
> 5 Memorial Dr, Building 5, Room 205.
> Bethesda, MD 20892. USA.
> Phone: 1-301-496-1592
> Fax: 1-301-496-9878
>
>
>
> On Mon, Mar 1, 2010 at 12:55 PM, Alec Wysoker <alecw@... 
> > wrote:
>> Release notes for Picard 1.15
>> 1 March 2010
>>
>> - Add code to write LiftOver Chain file.
>> - Add validation to LiftOver.
>> - In CommandLineParser, allow a special symbol to be used in  
>> parameter
>> names to force the option to be used even if coming from the params  
>> file
>> and would usually be overridden by the command line.
>> - Remove final from various public attribute-related SAMRecord  
>> methods.
>> I'm not sure why these were final to begin with.
>> - Added support for clipping to SamToFastq.
>> - Refactored cigar-clipping into a separate iterator for use  
>> outside of
>> command-line program CleanSam.
>> - Add utility to reconstruct reference bases from read with cigar  
>> and MD
>> string.  Not tested exhaustively, so please report any problems.
>> - Warn if paired read is not marked as first or second of pair.
>> - Add SamLocusIterator -- class for iterating over a SAM file by  
>> locus,
>> and getting pile of reads and offsets that align with locus.
>> - Add getter and setter for predicted median insert size to
>> SAMReadGroupRecord.
>> - Modified ValidateSamFile to return a non-zero code when  
>> validation fails.
>> - Change ValidateSamFile so that invalid value for @RG PI tag  
>> triggers
>> error message instead of exception.
>>
>> -Alec
>>
>> ------------------------------------------------------------------------------
>> Download Intel&#174; Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Samtools-help mailing list
>> Samtools-help@...
>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help



-- 
 The Wellcome Trust Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE.

Re: [Samtools-help] [Samtools-devel] Picard release 1.15

From: Keiran Raine <kr2@sa...> - 2010-03-02 08:32

Of course I have made an assumption that you used BWA here, so  
disregard if you have something that is mapping circular genomes  
correctly.

Keiran

On 2 Mar 2010, at 08:30, Keiran Raine wrote:

> Hi,
>
> What you are seeing is actually a bug/feature from BWA (which should
> not occur in the latest release).
>
> Your read is actually mapping to the last 33 bases of chrM and the
> first 3 bases of chr13_random.  Under the latest release this read
> would be marked as unmapped.
>
> Regards,
>
> Keiran Raine
> Senior Computer Biologist
> The Cancer Genome Project
> Ext: 2100
> kr2@...
>
> The Wellcome Trust Sanger Institute is operated by Genome Research
> Limited, a charity registered in England with number 1021457 and a
> company registered in England with number 2742969, whose registered
> office is 215 Euston Road, London, NW1 2BE.
>
>
>
> On 1 Mar 2010, at 19:38, Ivan Gregoretti wrote:
>
>> Hi again,
>>
>> This time I bring up an issue of bigger consequence that I just
>> spotted in this release. It has to do with boundary conditions on BAM
>> files. This is the example:
>>
>> I have a SAM record that goes like this
>>
>> zcat s_5.sam.gz | perl -wlne'if($_=~/HWI-EAS179_1:5:71:1726:6852*/)
>> {print}'
>> HWI-EAS179_1:5:71:1726:6852	0	chrM	16266	118	36M	*	0	0	
>> TAATCATACTCTATTACGCAATAAACATTAACAAGT	?>:;77/8<7A=BAABDDA=BB@...?
>> B<DB	XD:Z:34NN	SM:i:118
>>
>> As you can see, it is 36 bp record that starts at position 16266 and
>> it matches the forward strand of chrM.
>>
>> Now, the corresponding BAM file was successfully created like this
>>
>> samtools view -t mouse.mm9.genome -bS s_5.sam.gz > s_5.bam
>>
>> where an excerpt of 'cat mouse.mm9.genome' goes like this
>> ...
>> chrY	15902555
>> chrM	16299
>> chr13_random	400311
>> ...
>>
>> Notice that 36 + 16266 = 16302. Ugh.
>> ViewSAM.jar complains that this records is out of range and
>> terminates parsing.
>>
>> Now, mitochondrial DNA is linear!
>>
>> Ivan
>>
>>
>> Ivan Gregoretti, PhD
>> National Institute of Diabetes and Digestive and Kidney Diseases
>> National Institutes of Health
>> 5 Memorial Dr, Building 5, Room 205.
>> Bethesda, MD 20892. USA.
>> Phone: 1-301-496-1592
>> Fax: 1-301-496-9878
>>
>>
>>
>> On Mon, Mar 1, 2010 at 12:55 PM, Alec Wysoker <alecw@...
>>> wrote:
>>> Release notes for Picard 1.15
>>> 1 March 2010
>>>
>>> - Add code to write LiftOver Chain file.
>>> - Add validation to LiftOver.
>>> - In CommandLineParser, allow a special symbol to be used in
>>> parameter
>>> names to force the option to be used even if coming from the params
>>> file
>>> and would usually be overridden by the command line.
>>> - Remove final from various public attribute-related SAMRecord
>>> methods.
>>> I'm not sure why these were final to begin with.
>>> - Added support for clipping to SamToFastq.
>>> - Refactored cigar-clipping into a separate iterator for use
>>> outside of
>>> command-line program CleanSam.
>>> - Add utility to reconstruct reference bases from read with cigar
>>> and MD
>>> string.  Not tested exhaustively, so please report any problems.
>>> - Warn if paired read is not marked as first or second of pair.
>>> - Add SamLocusIterator -- class for iterating over a SAM file by
>>> locus,
>>> and getting pile of reads and offsets that align with locus.
>>> - Add getter and setter for predicted median insert size to
>>> SAMReadGroupRecord.
>>> - Modified ValidateSamFile to return a non-zero code when
>>> validation fails.
>>> - Change ValidateSamFile so that invalid value for @RG PI tag
>>> triggers
>>> error message instead of exception.
>>>
>>> -Alec
>>>
>>> ------------------------------------------------------------------------------
>>> Download Intel&#174; Parallel Studio Eval
>>> Try the new software tools for yourself. Speed compiling, find bugs
>>> proactively, and fine-tune applications for parallel performance.
>>> See why Intel Parallel Studio got high marks during beta.
>>> http://p.sf.net/sfu/intel-sw-dev
>>> _______________________________________________
>>> Samtools-help mailing list
>>> Samtools-help@...
>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>>
>>
>> ------------------------------------------------------------------------------
>> Download Intel&#174; Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Samtools-help mailing list
>> Samtools-help@...
>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>
>
>
> -- 
> The Wellcome Trust Sanger Institute is operated by Genome Research
> Limited, a charity registered in England with number 1021457 and a
> company registered in England with number 2742969, whose registered
> office is 215 Euston Road, London, NW1 2BE.
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-devel mailing list
> Samtools-devel@...
> https://lists.sourceforge.net/lists/listinfo/samtools-devel



-- 
 The Wellcome Trust Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE.

Re: [Samtools-help] [Samtools-devel] Picard release 1.15

From: Ivan Gregoretti <ivangreg@gm...> - 2010-03-02 14:44

Hello Keiran,

Indeed, I used Eland to map the sequence. Thank you for the insight.

Ivan

Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1016
Fax: 1-301-496-9878



On Tue, Mar 2, 2010 at 3:31 AM, Keiran Raine <kr2@...> wrote:
> Of course I have made an assumption that you used BWA here, so disregard if
> you have something that is mapping circular genomes correctly.
>
> Keiran
>
> On 2 Mar 2010, at 08:30, Keiran Raine wrote:
>
>> Hi,
>>
>> What you are seeing is actually a bug/feature from BWA (which should
>> not occur in the latest release).
>>
>> Your read is actually mapping to the last 33 bases of chrM and the
>> first 3 bases of chr13_random.  Under the latest release this read
>> would be marked as unmapped.
>>
>> Regards,
>>
>> Keiran Raine
>> Senior Computer Biologist
>> The Cancer Genome Project
>> Ext: 2100
>> kr2@...
>>
>> The Wellcome Trust Sanger Institute is operated by Genome Research
>> Limited, a charity registered in England with number 1021457 and a
>> company registered in England with number 2742969, whose registered
>> office is 215 Euston Road, London, NW1 2BE.
>>
>>
>>
>> On 1 Mar 2010, at 19:38, Ivan Gregoretti wrote:
>>
>>> Hi again,
>>>
>>> This time I bring up an issue of bigger consequence that I just
>>> spotted in this release. It has to do with boundary conditions on BAM
>>> files. This is the example:
>>>
>>> I have a SAM record that goes like this
>>>
>>> zcat s_5.sam.gz | perl -wlne'if($_=~/HWI-EAS179_1:5:71:1726:6852*/)
>>> {print}'
>>> HWI-EAS179_1:5:71:1726:6852     0       chrM    16266   118     36M     *
>>>       0       0
>>> TAATCATACTCTATTACGCAATAAACATTAACAAGT    ?>:;77/8<7A=BAABDDA=BB@...?
>>> B<DB    XD:Z:34NN       SM:i:118
>>>
>>> As you can see, it is 36 bp record that starts at position 16266 and
>>> it matches the forward strand of chrM.
>>>
>>> Now, the corresponding BAM file was successfully created like this
>>>
>>> samtools view -t mouse.mm9.genome -bS s_5.sam.gz > s_5.bam
>>>
>>> where an excerpt of 'cat mouse.mm9.genome' goes like this
>>> ...
>>> chrY    15902555
>>> chrM    16299
>>> chr13_random    400311
>>> ...
>>>
>>> Notice that 36 + 16266 = 16302. Ugh.
>>> ViewSAM.jar complains that this records is out of range and
>>> terminates parsing.
>>>
>>> Now, mitochondrial DNA is linear!
>>>
>>> Ivan
>>>
>>>
>>> Ivan Gregoretti, PhD
>>> National Institute of Diabetes and Digestive and Kidney Diseases
>>> National Institutes of Health
>>> 5 Memorial Dr, Building 5, Room 205.
>>> Bethesda, MD 20892. USA.
>>> Phone: 1-301-496-1592
>>> Fax: 1-301-496-9878
>>>
>>>
>>>
>>> On Mon, Mar 1, 2010 at 12:55 PM, Alec Wysoker <alecw@...
>>>>
>>>> wrote:
>>>> Release notes for Picard 1.15
>>>> 1 March 2010
>>>>
>>>> - Add code to write LiftOver Chain file.
>>>> - Add validation to LiftOver.
>>>> - In CommandLineParser, allow a special symbol to be used in
>>>> parameter
>>>> names to force the option to be used even if coming from the params
>>>> file
>>>> and would usually be overridden by the command line.
>>>> - Remove final from various public attribute-related SAMRecord
>>>> methods.
>>>> I'm not sure why these were final to begin with.
>>>> - Added support for clipping to SamToFastq.
>>>> - Refactored cigar-clipping into a separate iterator for use
>>>> outside of
>>>> command-line program CleanSam.
>>>> - Add utility to reconstruct reference bases from read with cigar
>>>> and MD
>>>> string.  Not tested exhaustively, so please report any problems.
>>>> - Warn if paired read is not marked as first or second of pair.
>>>> - Add SamLocusIterator -- class for iterating over a SAM file by
>>>> locus,
>>>> and getting pile of reads and offsets that align with locus.
>>>> - Add getter and setter for predicted median insert size to
>>>> SAMReadGroupRecord.
>>>> - Modified ValidateSamFile to return a non-zero code when
>>>> validation fails.
>>>> - Change ValidateSamFile so that invalid value for @RG PI tag
>>>> triggers
>>>> error message instead of exception.
>>>>
>>>> -Alec
>>>>
>>>>
>>>> ------------------------------------------------------------------------------
>>>> Download Intel&#174; Parallel Studio Eval
>>>> Try the new software tools for yourself. Speed compiling, find bugs
>>>> proactively, and fine-tune applications for parallel performance.
>>>> See why Intel Parallel Studio got high marks during beta.
>>>> http://p.sf.net/sfu/intel-sw-dev
>>>> _______________________________________________
>>>> Samtools-help mailing list
>>>> Samtools-help@...
>>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>>>
>>>
>>>
>>> ------------------------------------------------------------------------------
>>> Download Intel&#174; Parallel Studio Eval
>>> Try the new software tools for yourself. Speed compiling, find bugs
>>> proactively, and fine-tune applications for parallel performance.
>>> See why Intel Parallel Studio got high marks during beta.
>>> http://p.sf.net/sfu/intel-sw-dev
>>> _______________________________________________
>>> Samtools-help mailing list
>>> Samtools-help@...
>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>
>>
>>
>> --
>> The Wellcome Trust Sanger Institute is operated by Genome Research
>> Limited, a charity registered in England with number 1021457 and a
>> company registered in England with number 2742969, whose registered
>> office is 215 Euston Road, London, NW1 2BE.
>>
>>
>> ------------------------------------------------------------------------------
>> Download Intel&#174; Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Samtools-devel mailing list
>> Samtools-devel@...
>> https://lists.sourceforge.net/lists/listinfo/samtools-devel
>
>
>
> --
> The Wellcome Trust Sanger Institute is operated by Genome ResearchLimited, a
> charity registered in England with number 1021457 and acompany registered in
> England with number 2742969, whose registeredoffice is 215 Euston Road,
> London, NW1 2BE.
>

Re: [Samtools-help] Picard merging @PG issue

From: Alec Wysoker <alecw@br...> - 2010-03-01 22:35

Attachments: Message as HTML

Hi Folks,

We are assuming that the lack of objection to my proposal over a couple 
of weeks implies consent.  We're going to go ahead and implement these 
changes (optional PP and PN tags in PG header record) in Picard.  Can 
the spec be updated to reflect this?

Thanks, Alec

Goncalo Abecasis wrote:
> I very much like Alec's solution of chaining program tags with the PN and PP
> ids. It seems much better than using multiple tags in one line (harder to
> parse and digest) and also  straightforward to generate (a new program that
> modifies a BAM just has to point a new PN ID to the root of the previous
> sequence, which can otherwise remain unchanged).
>
> Goncalo
>
>
> -----Original Message-----
> From: Alec Wysoker [mailto:alecw@...] 
> Sent: Thursday, February 18, 2010 10:59 AM
> To: Richard Durbin
> Cc: samtools-help@...
> Subject: Re: [Samtools-help] Picard merging @PG issue
>
> Hi Folks,
>
> Has there been any response from Goncalo or Heng?  If no one objects, we
> would like to update the spec and implement this in Picard.
>
> Thanks, Alec
>
> Richard Durbin wrote:
>   
>> I quite like Alec's solution.  The complication is that when you apply 
>> a new program uniformly to a file, you will generate a new PG line for 
>> each current one (i.e. each one listed in a standard record).  But I 
>> think that is OK.  I would be in favour of requiring a PN, and using 
>> an arbitrary, e.g. normally int, identifier for the id.  I guess 
>> people want backward compatibility.  But if we get Picard, samtools 
>> and GATK implementing the PN and PP tags then that would go a long 
>> way.  I think I prefer this solution to listing multiple PG ids in 
>> each record, as Bob proposes, though I acknowledge that his proposal 
>> would work, so am open to continued consideration of it.  I'd quite 
>> like to know what Heng and Goncalo think.
>>
>> Richard
>>
>> Alec Wysoker wrote:
>>     
>>> Hi Folks,
>>>
>>> I think there is a way to avoid changing the basic syntax that will 
>>> accomplish the goal of recording an ordered list of programs that 
>>> have operated on alignments.  A new PP (previous program record) tag 
>>> could be added to the PG record.  This would contain the PG ID of the 
>>> previous program record.  A SAM record would store in its PG tag the 
>>> ID of the most recent PG record, but the chain of programs could be 
>>> discovered by following the chain of PP IDs.
>>>
>>> BTW, the proliferation of PG records that will eventually result 
>>> makes it even more important that we store program name separately so 
>>> that we can handle collisions of PG ID by arbitrarily changing the IDs.
>>>
>>> -Alec
>>>
>>> Bob Handsaker wrote:
>>>       
>>>> On 2/12/10 11:17 AM, Richard Durbin wrote:
>>>>  
>>>>         
>>>>> Let us think about the semantics of what we are trying to achieve.
>>>>>
>>>>> A sequence of commands has been applied to each record.  We want to 
>>>>> record that history.
>>>>>
>>>>> Clearly a set is not adequate to do that - we need something that 
>>>>> establishes the order.
>>>>>
>>>>> There seem to me two choices:
>>>>>      a) have one PG id per record, and allow there to be an ordered 
>>>>> sequence of multiple commands corresponding to that PG id
>>>>>      b) have multiple, ordered PG ids per record I prefer (a) over 
>>>>> (b) because it keeps the per-record information
>>>>> simple: there are millions of records and we want them to be 
>>>>> compact and simple to parse and operate on.
>>>>>        
>>>>>           
>>>> It seems to me that (b) is not necessarily a bad idea.
>>>> We could extend the syntax of the PG tag on the alignment to record 
>>>> allow multiple IDs, for example separated by a semi-colon:
>>>>      PG:bwa;gatk
>>>> Since in many cases the number of tag values in the file would be a 
>>>> small set, this might not have much overhead after bam compression.
>>>> Or if it was a problem in practice, it seems simple to design a more 
>>>> efficient and backwards compatible internal representation in bam.
>>>> This solution has the advantage that it is very simple and clear 
>>>> what to do if you are reprocessing existing records.
>>>> -Bob
>>>>  
>>>>         
>>>>> Given (a) we have two options again
>>>>>      (a1) have one header line per PG id, with the sequence of 
>>>>> multiple commands on that line
>>>>>      (a2) have multiple header lines per PG id, with the order of 
>>>>> those lines meaningful Of these I prefer (a1) over (a2) for 
>>>>> multiple reasons: first it means that PG ids are unique on header 
>>>>> lines, which seems natural and expected.  Second, even if there 
>>>>> were multiple header lines, I don't very much like the semantics 
>>>>> that the order in the header of those lines is meaningful - I think 
>>>>> that header semantics are that header lines are a set.
>>>>> On the other hand, generic header code would be expected to move 
>>>>> header lines around wholesale, not reorder within them (merge code 
>>>>> could remap the ID if the following text was not identical).  
>>>>> Whatever code interprets the PG header lines can know that within 
>>>>> these the order of items is meaningful.
>>>>>
>>>>> Having established (a1), I agree that flat subfields don't look so 
>>>>> nice.  I would prefer some syntax like:
>>>>>
>>>>> PG<pgid>    '('<name>  [VN:version] [CL:arguments]')'+
>>>>>
>>>>> i.e. incorporate the information about each command in parentheses 
>>>>> (or some other sort of paired brackets such as {} or [].
>>>>>
>>>>> Richard
>>>>>
>>>>> John Marshall wrote:
>>>>>       
>>>>>           
>>>>>> On 12 Feb 2010, at 13:52, Richard Durbin wrote:
>>>>>>           
>>>>>>             
>>>>>>> It looks like we are getting convergence.  As I said in a message 
>>>>>>> that I think crossed this, it seems to me that what we want is 
>>>>>>> for each record the history of commands that were applied to make 
>>>>>>> it, and modify it.
>>>>>>>                
>>>>>>>               
>>>>>> [namely...]
>>>>>>           
>>>>>>             
>>>>>>>> PG  ID:pgid1 CM:bwa  VN:0.5.5  CM:GATK VN:v2.2.16 
>>>>>>>> CL:Covariates=[ReadGroupCovariate,...
>>>>>>>>                    
>>>>>>>>                 
>>>>>>           
>>>>>>             
>>>>>>> I think that each each record should have one PG id, and each PG 
>>>>>>> id should allow a sequence of commands.
>>>>>>> Is that what you are proposing below?  Could there be multiple 
>>>>>>> PN, VN, CL entries on one PG header line?
>>>>>>>                
>>>>>>>               
>>>>>> While the spec does not explicitly forbid header fields appearing 
>>>>>> more than once in the same header, some current implementations 
>>>>>> do: it will fail Picard validation.  IMHO Picard has the right 
>>>>>> idea and the spec should be tightened to forbid this, as it 
>>>>>> already does for optional alignment auxiliary fields.
>>>>>>
>>>>>> At the moment there is no meaning to the order in which header 
>>>>>> fields appear, and tools might well ignore ordering when doing 
>>>>>> things like determining whether two headers (e.g., from different 
>>>>>> files to be
>>>>>> merged) are equivalent.  In fact, a header's fields currently form 
>>>>>> a mathematical set, which has simple well-understood properties.  
>>>>>> Under this proposal, they would be much more complicated: a 
>>>>>> multiset with at least some meaningful ordering (or else how would 
>>>>>> you know which VN went with which CM?).  We'd either have to 
>>>>>> formulate complicated rules for when headers are equivalent etc, 
>>>>>> or the tools would be unable to unify differently presented 
>>>>>> headers that a human could see were clearly equivalent.
>>>>>>
>>>>>> Both of those sound much worse than the current simplicity.
>>>>>>
>>>>>>      John
>>>>>>            
>>>>>>             
>>>>>        
>>>>>           
>>>> --------------------------------------------------------------------
>>>> ----------
>>>>
>>>> SOLARIS 10 is the OS for Data Centers - provides features such as 
>>>> DTrace, Predictive Self Healing and Award Winning ZFS. Get Solaris 
>>>> 10 NOW http://p.sf.net/sfu/solaris-dev2dev
>>>> _______________________________________________
>>>> Samtools-help mailing list
>>>> Samtools-help@...
>>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>>>   
>>>>         
>>> ---------------------------------------------------------------------
>>> ---
>>>
>>> ---------------------------------------------------------------------
>>> ---------
>>>
>>> SOLARIS 10 is the OS for Data Centers - provides features such as 
>>> DTrace, Predictive Self Healing and Award Winning ZFS. Get Solaris 10 
>>> NOW http://p.sf.net/sfu/solaris-dev2dev
>>> ---------------------------------------------------------------------
>>> ---
>>>
>>> _______________________________________________
>>> Samtools-help mailing list
>>> Samtools-help@...
>>> https://lists.sourceforge.net/lists/listinfo/samtools-help
>>>   
>>>       
>>
>>     
>
> ----------------------------------------------------------------------------
> --
> Download Intel&#174; Parallel Studio Eval Try the new software tools for
> yourself. Speed compiling, find bugs proactively, and fine-tune applications
> for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Samtools-help mailing list
> Samtools-help@...
> https://lists.sourceforge.net/lists/listinfo/samtools-help
>
>
>
>

Re: [Samtools-help] SAM format and strandedness

From: Angie Hinrichs <angie@so...> - 2010-03-01 18:17

Attachments: Message as HTML

> flag >= 16 -> strand - 

For single-ended reads that is right, but I hope nobody gets the idea to try that with paired-end reads. Also, I have never seen the topmost flag bits (BAM_FSECONDARY, BAM_FQCFAIL, BAM_FDUP) used in practice, but they would break that rule for single-ended reads. I find hexadecimal to be a much more clear representation of bitwise flags... there is only one hex digit to examine for a given bit position, and I never have to do math with numbers greater than f. :) 

Angie 

----- "Davide Cittaro" <davide.cittaro@...> wrote: 
> From: "Davide Cittaro" <davide.cittaro@...> 
> To: "samtools-help" <samtools-help@...> 
> Sent: Sunday, February 28, 2010 12:40:33 AM GMT -08:00 US/Canada Pacific 
> Subject: Re: [Samtools-help] SAM format and strandedness 
> 
> 
> 

On Feb 27, 2010, at 12:24 PM, Davide Cittaro wrote: 

> flag < 16 -> strand + 
> flag >= 16 -> strand - 
> 
> (this is if one can't check bitwise and with 0x04) 
> 
> 

> 
I meant 0x0010 obviously... 

> 
d 
> 

> 
> ------------------------------------------------------------------------------ Download Intel® Parallel Studio Eval Try the new software tools for yourself. Speed compiling, find bugs proactively, and fine-tune applications for parallel performance. See why Intel Parallel Studio got high marks during beta. http://p.sf.net/sfu/intel-sw-dev 
> _______________________________________________ Samtools-help mailing list Samtools-help@... https://lists.sourceforge.net/lists/listinfo/samtools-help

[Samtools-help] samtools' pileup format

From: Mehta, Perdeep <Perdeep.Mehta@ST...> - 2010-03-01 16:16

Attachments: Message as HTML

Hi,

I have a file generated by samtools pileup option.  The output is fine and everything looks ok, except at few places annotation does not make sense or I may be missing something.  According to samtools doc for pileup, a symbol '^' marks the start of a read segment which is a contiguous subsequence on the read separated by 'N/S/H' CIGAR operations. A symbol '$' marks the end of a read segment.

A preceding read in the following example ends at 590115 with '$' mark, but I don't see a '^' at the start of next read beginning at 631553, but it does end at 631582 with '$' mark.  It does not seems to be consistent if there is single or multiple mapped reads to a location. If formatting of the example gets screwed up, I have attached a word doc.

Any clues or thoughts will help.

Thanks,
perdeep

-
-
chr10   590112  G       G       84      0       51      19      ,,,,,,,,,,,,,,,,,,, 7:@.?<:79;;9(4614<7
chr10   590113  T       T       84      0       51      19      ,,,,,,,,,,,,,,,,,,,     6941:@:<0.+0;71/)B2
chr10   590114  G       G       84      0       51      19      ,$,$,$,$,,$,$,$,$,$,$,$,$,$,,$,$,,$     A3815501;6;5181;1<=
chr10   590115  G       G       18      0       54      3       ,$,$,$  )((
chr10   631553  a       A       4       0       1       1       .       2
chr10   631554  a       A       4       0       1       1       .       9
chr10   631555  a       G       4       4       1       1       G       2
chr10   631556  a       A       4       0       1       1       .       2
chr10   631557  a       A       4       0       1       1       .       =
chr10   631558  a       G       4       4       1       1       G       2
chr10   631559  a       A       4       0       1       1       .       1
chr10   631560  a       A       4       0       1       1       .       @
chr10   631561  a       A       4       0       1       1       .       2
chr10   631562  a       A       4       0       1       1       .       @
chr10   631563  a       A       4       0       1       1       .       1
chr10   631564  a       A       4       0       1       1       .       9
chr10   631565  a       A       4       0       1       1       .       @
chr10   631566  a       A       4       0       1       1       .       ;
chr10   631567  a       A       4       0       1       1       .       (
chr10   631568  a       A       4       0       1       1       .       :
chr10   631569  a       A       4       0       1       1       .       @
chr10   631570  T       T       4       0       1       1       .       2
chr10   631571  T       T       4       0       1       1       .       >
chr10   631572  A       A       4       0       1       1       .       +
chr10   631573  G       G       4       0       1       1       .       (
chr10   631574  A       A       4       0       1       1       .       >
chr10   631575  T       T       4       0       1       1       .       5
chr10   631576  G       T       4       4       1       1       T       >
chr10   631577  T       T       4       0       1       1       .       0
chr10   631578  A       A       4       0       1       1       .       5
chr10   631579  C       C       4       0       1       1       .       0
chr10   631580  T       T       4       0       1       1       .       0
chr10   631581  C       C       4       0       1       1       .       8
chr10   631582  A       A       4       0       1       1       .$      8
chr10   631914  a       A       12      0       12      1       ^-,     :
chr10   631915  a       A       12      0       12      1       ,       :
chr10   631916  a       A       12      0       12      1       ,       @
chr10   631917  a       A       12      0       12      1       ,       @
-
-



  ________________________________
Email Disclaimer: http://www.stjude.org/emaildisclaimer

[Samtools-help] export2sam doubts...

From: Gopu <gopu_36@ya...> - 2010-03-01 15:42

Attachments: Message as HTML

Hi,

I have few queries about samtools.

1. I am using eland mapping output and start using export2sam.pl. All the PF reads from export are being used for down stream analysis. 
2. How the uniquely mapped and multiple mapped are handled during pileup command? Please confirm whether only unique hits from export.txt are used for pileUp command or all PF reads.
3. The extended CIGAR column always shows 35M (ie the length of the read). How did the mismatch information would be incorporated? The column 15 of export contains the match descriptor information.

Thanks. 



      New Email names for you! 
Get the Email name you&#39;ve always wanted on the new @ymail and @rocketmail. 
Hurry before someone else does!
http://mail.promotions.yahoo.com/newdomains/sg/

[Samtools-help] Troubleshoot samtools in Windows environment

From: Li Zhihui <zhihui@ge...> - 2010-02-26 10:10

Attachments: Message as HTML Message as HTML

Next Messages