Question: Are secondary alignments (SAM) produced by STAR aligner linked with MAPQ score?
0
gravatar for Claudio
2.4 years ago by
Claudio0
Claudio0 wrote:

Hi everyone,

I am currently working with paired-end RNA-Seq data and I have problems in understanding some details of the SAM content. Let us use for instance one BAM file available in ENCODE

EXPERIMENT: https://www.encodeproject.org/experiments/ENCSR000AEW/

DIRECT DOWNLOAD: https://www.encodeproject.org/files/ENCFF602BYA/@@download/ENCFF602BYA.bam

The BAM file has been produced using STAR. From the header it is possible to retrieve the executed command:

@CO user command line: STAR --genomeDir out --readFilesIn ENCFF001ROL.fastq.gz ENCFF001ROK.fastq.gz --readFilesCommand zcat --runThreadN 8 --genomeLoad NoSharedMemory --outFilterMultimapNmax 20 --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --outFilterMismatchNmax 999 --outFilterMismatchNoverReadLmax 0.04 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 --outSAMheaderCommentFile COfile.txt --outSAMheaderHD @HD VN:1.4 SO:coordinate --outSAMunmapped Within --outFilterType BySJout --outSAMattributes NH HI AS NM MD --outSAMtype BAM SortedByCoordinate --quantMode TranscriptomeSAM --sjdbScore 1 --limitBAMsortRAM 60000000000

From the ENCFF602BYA.bam file I filtered all the reads having read name "DFDF8JF1:270:C19DCACXX:2:1216:1336:52997". It produces 6 entries.

DFDF8JF1:270:C19DCACXX:2:1216:1336:52997    99  chr11   83874510    1   45M87726N54M65533N2M    =   84148457    274048  CAGTGTAATTTCTTCAAATTCATATTCAATTTCTGTCCCATTGACATAAGGAATTGTGTCCAAAGTATCTGTGTTGACAATTATAGGAGCAGGACTGGCCT   @@@FFDDFHHDFHGHGGJGIJJCIICII@FFHGIEHHEGFFG<@F<<FHDHGFHDGBFGHG@FEG??EFEHFGBC=)@@77BGGI4.?#############   NH:i:3  HI:i:1  AS:i:195    NM:i:0  MD:Z:101

DFDF8JF1:270:C19DCACXX:2:1216:1336:52997    355 chr11   83874510    1   45M87726N54M186096N2M   =   84148457    274048  CAGTGTAATTTCTTCAAATTCATATTCAATTTCTGTCCCATTGACATAAGGAATTGTGTCCAAAGTATCTGTGTTGACAATTATAGGAGCAGGACTGGCCT   @@@FFDDFHHDFHGHGGJGIJJCIICII@FFHGIEHHEGFFG<@F<<FHDHGFHDGBFGHG@FEG??EFEHFGBC=)@@77BGGI4.?#############   NH:i:3  HI:i:2  AS:i:195    NM:i:0  MD:Z:101

DFDF8JF1:270:C19DCACXX:2:1216:1336:52997    355 chr11   83874510    1   45M87726N56M    =   84148457    274048  CAGTGTAATTTCTTCAAATTCATATTCAATTTCTGTCCCATTGACATAAGGAATTGTGTCCAAAGTATCTGTGTTGACAATTATAGGAGCAGGACTGGCCT   @@@FFDDFHHDFHGHGGJGIJJCIICII@FFHGIEHHEGFFG<@F<<FHDHGFHDGBFGHG@FEG??EFEHFGBC=)@@77BGGI4.?#############   NH:i:3  HI:i:3  AS:i:194    NM:i:0  MD:Z:101

DFDF8JF1:270:C19DCACXX:2:1216:1336:52997    147 chr11   84148457    1   101M    =   83874510    -274048 TGTCCACTAACCTGCGACCCAGCTTCTTAGCATACCAGATAGAGGCAAACATAGCGATGCCAATGGAAGGAATTCACAATGAAGCTGCCGCACAATCAACT   ###########BB;3(=BCC@>3>;):?7.)?>A;CAGIHDC=4.)=.)>A=AFGDD90*?4FBAGGGEC?:**9?9*H9EA2+2)E@8C:=ADD?DD??@   NH:i:3  HI:i:1  AS:i:195    NM:i:2  MD:Z:3G3G93

DFDF8JF1:270:C19DCACXX:2:1216:1336:52997    403 chr11   84148457    1   101M    =   83874510    -274048 TGTCCACTAACCTGCGACCCAGCTTCTTAGCATACCAGATAGAGGCAAACATAGCGATGCCAATGGAAGGAATTCACAATGAAGCTGCCGCACAATCAACT   ###########BB;3(=BCC@>3>;):?7.)?>A;CAGIHDC=4.)=.)>A=AFGDD90*?4FBAGGGEC?:**9?9*H9EA2+2)E@8C:=ADD?DD??@   NH:i:3  HI:i:2  AS:i:195    NM:i:2  MD:Z:3G3G93

DFDF8JF1:270:C19DCACXX:2:1216:1336:52997    403 chr11   84148457    1   101M    =   83874510    -274048 TGTCCACTAACCTGCGACCCAGCTTCTTAGCATACCAGATAGAGGCAAACATAGCGATGCCAATGGAAGGAATTCACAATGAAGCTGCCGCACAATCAACT   ###########BB;3(=BCC@>3>;):?7.)?>A;CAGIHDC=4.)=.)>A=AFGDD90*?4FBAGGGEC?:**9?9*H9EA2+2)E@8C:=ADD?DD??@   NH:i:3  HI:i:3  AS:i:194    NM:i:2  MD:Z:3G3G93

I undestand most of the information reported here. Because it is paired-end data, each read name is present at least twice. Once per mate. From this example there is one primary alignment (first/fourth entry) and two secondary alignments (second/fifth entry, third/sixth entry). The multiple alignment is expressed by the different CIGAR score in the first three entries.

I am investigating where reads split across the genome, so to identify those exons that are expressed in the same isoform (de novo analysis). Now my questions are three:

  1. Out of the three multiple alignments, should I take the primary one while considering the coverage?
  2. Is the MAPQ score related with the number of multiple alignments? In the STAR documentation, "1" means "Maps to 4-9 locations in the target", but here I see only three multiple alignment. Does MAPQ involve also the CIGAR information?
  3. "AS" SAM tag express the "Alignment score generated by aligner", how can I put the scores provided in this example in context?

Thanks in advance for the help

Claudio

rna-seq alignment forum • 1.2k views
ADD COMMENTlink modified 2.4 years ago • written 2.4 years ago by Claudio0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1357 users visited in the last hour