Hello!

I am working on RNA Seq. I mapped the reads to a genome using bwa. Now I
want to get the count tables (number of reads per genomic feature), and one
of the steps include conversion of .sam files into sorted .bam files, and
creating a GRangesList out of sorted .bam files. Before the last mentioned
step, I also applied readGappedAlignments to read in the sorted .bam files.
Code, tail of tested .sam file, and contest of variable bwa_output are
following:

*Code*:

fcatusDb=makeTranscriptDbFromBiomart(biomart="ensembl",
dataset="fcatus_gene_ensembl")
tx=transcriptsBy(fcatusDb)

bwa_output=readGappedAlignments(file_name1, format="BAM")

bwa_grglist=grglist(bwa_output)

rc=countOverlaps(tx, bwa_grglist[countOverlaps(bwa_grglist,tx)==1])



*test.sam*:

@SQ     SN:GeneScaffold_295     LN:1999586
@SQ     SN:GeneScaffold_5264    LN:2008247
@SQ     SN:GeneScaffold_3594    LN:2061208
@SQ     SN:GeneScaffold_3107    LN:2081283
HWUSI-EAS475:1:1:9547:3828#0    0       GeneScaffold_1986       3069    0
    41M     *       0       0
CCACATCTCTCCCAGTTTCTTTGCAACATCACCAATGGATA       *       XT:A:R  NM:i:0
 X0:i:3  X1:i:1  XM:i:0  XO:i:0  XG:i:0  MD:Z:41
XA:Z:scaffold_214963,-224762,41M,0;GeneScaffold_4609,+41298,41M,0;GeneScaffold_886,-415265,41M,1;
HWUSI-EAS475:1:1:5401:3825#0    16      scaffold_96256  459     0       41M
    *       0       0       GTTGGCTCTCAGGGGTACATGATGAACTGGGGAAGGAGTAT
*       XT:A:R  NM:i:0  X0:i:3      X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:41
XA:Z:scaffold_166015,+60068,41M,0;scaffold_5553,+4693,41M,0;
HWUSI-EAS475:1:1:8428:7999#0    0       GeneScaffold_4103       98907   37
     41M     *       0       0
CCTGCACGTTCATTGTGTGTGTCTTGAGTTGATTTGTCACC       *       XT:A:U  NM:i:1
 X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:11T29
HWUSI-EAS475:1:1:5814:7999#0    0       scaffold_75884  112     14      41M
    *       0       0       CCCAGAAGTCCCAGGTGTTCCCATTCTGTGATCAGCACAGA
*       XT:A:U  NM:i:0  X0:i:1      X1:i:8  XM:i:0  XO:i:0  XG:i:0  MD:Z:41
XA:Z:scaffold_99844,+4698,41M,1;scaffold_20760,+365,41M,1;scaffold_95772,+1923,41M,1;scaffold_16284,-3782,41M,1;scaffold_22266,-52,41M,1;scaffold_82279,-5356,41M,1;scaffold_75127,-757,41M,1;scaffold_82292,-3298,41M,1;



*bwa_output:*
*
*
*
GappedAlignments with 4 alignments and 0 elementMetadata values:
               seqnames strand       cigar    qwidth     start       end
                  <Rle>  <Rle> <character> <integer> <integer> <integer>
  [1]    scaffold_75884      +         41M        41       112       152
  [2]    scaffold_96256      -         41M        41       459       499
  [3] GeneScaffold_1986      +         41M        41      3069      3109
  [4] GeneScaffold_4103      +         41M        41     98907     98947
          width      ngap
      <integer> <integer>
  [1]        41         0
  [2]        41         0
  [3]        41         0
  [4]        41         0
  ---


*
*
*

According the manual ("Multi-reads won't receive any special treatment i.e.
the various SAM/BAM records corresponding to a multi-read will show up in
the GappedAlignments object as if they were coming from different/unrelated
queries.") and taken into account previously mentioned, I would be grateful
if someone can confirm me that readGappedAlignments (output of bwa gives
.sam files which contain all information about mapping locations of one
read in one line) will only take the first hit, no matter what is following
(whether that read is mapped to several equally scored positions...).

If this is the case, would you filter the .sam files for uniquely mapped
reads before reading it via readGappedAlignments? Or will you let the
countOverlaps do it?



Thank you!

	[[alternative HTML version deleted]]

