[BioC] Rsamtools question

[Valerie Obenchain]

On 10/04/2011 08:58 AM, Valerie Obenchain wrote:
> Hi Alpesh,
>
> One approach would be to
>
> I. Read in snps with readGappedAlignments :
>
> Convert your sam file to a bam using samtools from the command line,
>
>     samtools view -bS -o  yourfile.bam  yourfile.sam

or from within R -- Rsamtools::asBam

>
> and read  into R with readGappedAlignments
>
>     library(GenomicFeatrures)
>     ?readGappedAlignments
>
>
> II. Get gene ID :
>
> Use an org annotation package to get the gene id you are interested 
> in. There are many ID types available (Entrez, Ensembl, RefGene, 
> etc.). You will use this ID to identify the gene of interest in the 
> TxDb package (below).
> Many org packages are available, see all starting with "org",
>
>     http://bioconductor.org/packages/devel/data/annotation/
>
> For the human package,
>
>     library(org.Hs.eg.db)
>     ?org.Hs.eg.db
>
>
> III. Get the genomic coordinates for the gene of interest :
> The TxDb packages have genomic coordinates as well as transcript, exon 
> and gene ID's. See packages starting with "TxDb",
>
>     http://bioconductor.org/packages/devel/data/annotation/
>
> For example, the human hg19 known genes from UCSC are in the following 
> package. Since this is a UCSC package the geneID's will be RefGene ID's.
>
>     library(GenomicFeatures)
>     library(TxDb.Hsapiens.UCSC.hg19.knownGene)
>     ls(2)  # to see the name of the txdb object
>
> If you want to make your own txdb see
>
>     ?makeTranscriptDb
>
> There are several functions available for organizing you data by 
> exons, transcripts or gene, see
>
>     ?transcriptsBy
>
>
> IV. Counting
>
> You now have the genomic coordinates of your gene (from txdb package) 
> and your snps (read in with readGappedAlignments).
> Several counting options are available,
>
>     library(GenomicRanges)
>     ?summarizeOverlaps # avoids double counting, available in R-devel 
> version
>     ?countOverlaps
>     ?findOverlaps

also Rsamtools::applyPileups which gives bases over genomic coordinates

Valerie
>
>
> Valerie
>
>
>
> On 10/02/2011 08:28 AM, Alpesh Querer wrote:
>> Hi All,
>>
>> I am a newbie with R and Bioconductor. I have .sam files and I want 
>> to count
>> the number of reads mapping to particular SNP locations on the .fa gene.
>> Also I want to avoid double counting the SNPs which are covered by 
>> the same
>> 36bp read.
>>   Any help is highly appreciated.
>>
>> Thanks,
>> Alpesh
>>
>>     [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives: 
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: 
> http://news.gmane.org/gmane.science.biology.informatics.conductor