Question: How to find the unique list of mismatches of reads with the reference genome
0
gravatar for Sama
11 days ago by
Sama30
Sama30 wrote:

Hi,

I have Illumina reads and a reference genome. I aligned reads to the genome and I have .bam and .sam files. I want the list of the mismatches of the reads with the reference genome but if one mismatch happened in multiple reads, I want it to be reported only once. For example if two reads have the mismatch of A->G on position 1000 of the reference genome, I want something like this: G,1000,2

I want the output file to be something like this:

Mismatch,Position,Frequency

G,1000,2

A,3278,7

G,78732,5

T,89783,8

C,87494,8

A,13732,5

...

Is there any tool available for doing that?

alignment next-gen genome • 122 views
ADD COMMENTlink modified 11 days ago by pfs90 • written 11 days ago by Sama30
3
gravatar for pfs
11 days ago by
pfs90
USA/Boston
pfs90 wrote:

It looks like you need to implement a variant calling step and then do some data parsing to get your desired out put format. See http://www.htslib.org/workflow/

ADD COMMENTlink written 11 days ago by pfs90
1
gravatar for Pierre Lindenbaum
11 days ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum99k wrote:

using sam2tsv: http://lindenb.github.io/jvarkit/Sam2Tsv.html


   $ java -jar dist/sam2tsv.jar -R ref.fa in.bam | awk '$8!=$5 && ($9=="M" || $9=="X")' | cut -f 3,4,5 | sort | uniq -c


      2 ref2    0   A
      1 ref2    0   C
      1 ref2    0   G
      2 ref2    0   T
      3 ref2    10  A
      2 ref2    10  T
      3 ref2    11  A
      1 ref2    11  C
      1 ref2    11  G
      1 ref2    12  A
      2 ref2    12  C
      2 ref2    12  T
      4 ref2    13  A
      2 ref2    13  C
      3 ref2    14  A
      1 ref2    14  C
      1 ref2    14  G
      1 ref2    14  T
      5 ref2    15  A
      1 ref2    15  T
      1 ref2    16  A
      1 ref2    16  C
      2 ref2    16  G
      2 ref2    16  T
      4 ref2    17  A
      1 ref2    17  C
      1 ref2    17  T
      4 ref2    18  A
      2 ref2    18  G
      3 ref2    19  A
      1 ref2    19  C
      1 ref2    19  G
      1 ref2    19  T
ADD COMMENTlink written 11 days ago by Pierre Lindenbaum99k

Thank you for your reply. I have a question. What are first and third columns?

ADD REPLYlink written 10 days ago by Sama30
1

in the final output ? occurence/chromosome/position/alt

ADD REPLYlink written 10 days ago by Pierre Lindenbaum99k

This is 5 first lines of running your command on my data:

3242056 Escherichia 0 A

3435134 Escherichia 0 C

3169976 Escherichia 0 G

85 Escherichia 0 N

2952800 Escherichia 0 T

It is weird if the frequency of A in first position is 3242056.

Also the third column numbers range is 0-100 but the reference genome is much longer than 100.

ADD REPLYlink modified 10 days ago • written 10 days ago by Sama30

Also the third column numbers range is 0-100 but the reference genome is much longer than 100.

I wonder if the offsets of the columns in my 'cut' are the correct one, may be I've used the position in the read. Can you please check (i'm away from my tools) with

 java -jar dist/sam2tsv.jar -R ref.fa in.bam | more
ADD REPLYlink written 10 days ago by Pierre Lindenbaum99k
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: 1417 users visited in the last hour