SEQanswers (
-   Bioinformatics (
-   -   SHRiMP 2 output problem (

rahilsethi 10-09-2012 10:08 AM

SHRiMP 2 output problem

I used SHRiMP 2.2.2 for mapping (gapped mapping) SOLiD Paired End data to hg19 genome reference. The SAM output file does not have any value at QV column, instead has "*" as a placeholder. For example:

271_822_1591_F 163 chrM 1 255 35M = 146 194 GATCACAGGTCTATCACCCTATTAACCACTCACGG * AS:i:310 NM:i:0 CS:Z:T12321112012233211002330301311222130 CM:i:2 XX:Z:GATCACAGGTCTATCACCCTATTAACcACTCaCGG
86_1077_1408_F 163 chrM 1 255 3H32M = 126 175 GATCACAGGTCTATCACCCTATTAACCACTCA * AS:i:320 NM:i:0 CS:Z:T33102321112012233211002330301011221 CM:i:0 XX:Z:GATCACAGGTCTATCACCCTATTAACCACTCA

I understand this because SHRiMP 2 never asks for .qual filename/s when put on run. Now almost all the variant calling algorithms such as SNVer, VarScan, GATK etc. require base/color quality values for evaluation of their variant call and end up giving error. Even the mpileup/bcftools seem to be giving abnormal results with almost all of the variant calls as In-Dels. Is there anyway I can make SHRiMP2 give out value to QV column in output so that I can run the variant calling tools on it?

My SHRiMP 2.2.2 run was:

gmapper-cs \
--max-alignments 800 \
-N 12 \
--pair-mode opp-in \
--isize 0,250 \
--all-contigs \
--single-best-mapping \
--local \
--sam \
-1 trimmed_read_F3_50bp.csfasta -2 read_F5.csfasta \
hg19.ucsc.fa > output_trimmed_50bp.sam 2> output_trimmed_50bp.log

Thanks in advance

nilshomer 10-09-2012 03:06 PM

The *csfasta files do not have any quality information, you probably need to input the *qual files as well.

rahilsethi 10-10-2012 05:05 AM

Hi nilshomer,

Thanks for the suggestion. Although SHRiMP, in general, does not give any separate option for .qual files in the run, I did try to incorporate qual files along with csfasta; but it prompts and error, indicating that it does not recognize the .qual format.
One suggestion has been recommended to me by SNVer authors is to include a dummy value in QUAL field of SHRiMP SAM output such as "IIIII" (size of SEQ length) for all mapped reads so as to generate variant calls without runtime error. But then, that would affect the filter process of the variants.
I have also come across the variant calling tool VARiD that shows an example of working on SHRiMP mapping output. I might give a try to VARiD. Do you know how good is VARiD tool for variant calling?

All times are GMT -8. The time now is 04:18 AM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.