Heads up:
We’re moving the GATK website, docs and forum to a new platform. Read the full story and breakdown of key changes on this blog.
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra

Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

Duplicate columns in ReadLengthDistribution

lindenblindenb FranceMember ✭✭
edited September 2016 in Ask the GATK team

Hi GATK team,
I've run the following command:

java -jar gatk/3.6.0/GenomeAnalysisTK.jar -T ReadLengthDistribution -R /commun/data/pubdb/broadinstitute.org/bundle/1.5/b37/human_g1k_v37.fasta -L "2:237200000-237210000" -o out.txt  -I  bam.list -rf NotPrimaryAlignment -rf UnmappedRead -rf FailsVendorQualityCheck -rf DuplicateRead

in out.txt: the header starts with

#:GATKTable:ReadLengthDistribution:Table of read length distributions
readLength  10H0268  12H0162  12H0162  12H0300  12H0300  12H0356  10H0192  10H0027  10H0027  11H0493  11H0493  14H0292  13H0199  12H0447  12H0447  12H0031  12H0031  11H0404  11H0404  10H0227  12H0473  10H0469  14H0292  10H0469  13H0490  13H0490  13H0025  13H0025  12H0341  12H0341  12H0098  12H0098  12H0604  12H0604

as you can see, there are some duplicate columns e.g:12H0300 (wich don't always have the same value in the table.

my bam.list contains only 20 bam:

$ wc  bam.list 
  20   20 2520 bam.list

looking at the BAM headers, I can see that there are some samples that have more than one read-grou

$ cat  bam.list  | while read F; do echo "###"; samtools view -H "${F}" | grep "^@RG" ; done 
@RG ID:p20  PL:ILLUMINA PU:1    LB:13H0199  SM:13H0199
@RG ID:p12  PL:ILLUMINA PU:1    LB:12H0300  SM:12H0300
@RG ID:p13  PL:ILLUMINA PU:1    LB:12H0300  SM:12H0300
@RG ID:p21  PL:ILLUMINA PU:1    LB:12H0447  SM:12H0447
@RG ID:p22  PL:ILLUMINA PU:1    LB:12H0447  SM:12H0447

I suppose ReadLengthDistribution group the data by group ID and display the SM tag. Am I right ? Do you consider this as a bug or as a feature ? :-)
Is there a way to group the data by SM ?


Best Answer


Sign In or Register to comment.