search

walaj / VariantBam

Filtering and profiling of next-generational sequencing data using region-specific rules
Other
74 stars 10 forks source link

Soft-clipping with linked-region filter #17

Open daynefiler opened 5 years ago

daynefiler commented 5 years ago

Hi, first -- thank you for your software! I am working with paired end sequencing on a pool of fragments that contain a large proportion of fragments shorter than the reads. This leads to soft-clipping in the (proper) alignment. Unfortunately, when filtering paired reads using the -l (linked region) flag, proper read pairs get included that do not overlap the given locus.

I've provided a MWE:

## github.sam ##
@HD VN:1.5  SO:coordinate
@SQ SN:NC_000001.11 LN:248956422
@RG ID:S3_L002  SM:S3   PL:Illumina CN:UNC
@RG ID:S3_L003  SM:S3   PL:Illumina CN:UNC
@PG ID:bwa  PN:bwa  VN:0.7.15-r1140 CL:bwa mem -t 8 -R @RG\tID:S3_L002\tSM:S3\tPL:Illumina\tCN:UNC /projects/sequence_analysis/vol4/bpow/references/GCF_000001405.26_GRCh38_genomic.fna inputs/S3/L002/S3_L002_R2.fastq.gz inputs/S3/L002/S3_L002_R1.fastq.gz
@PG ID:MarkDuplicates   VN:2.9.2-SNAPSHOT   CL:picard.sam.markduplicates.MarkDuplicates INPUT=[mapped/S3.L002.bam] OUTPUT=markdup/S3.L002.bam METRICS_FILE=markdup/S3.L002.metrics.txt REMOVE_DUPLICATES=false TMP_DIR=[/scratch/dfiler] CREATE_INDEX=true    MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 TAG_DUPLICATE_SET_MEMBERS=false REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag ASSUME_SORTED=false DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json    PN:MarkDuplicates
@PG ID:bwa-4A25094  PN:bwa  VN:0.7.15-r1140 CL:bwa mem -t 8 -R @RG\tID:S3_L003\tSM:S3\tPL:Illumina\tCN:UNC /projects/sequence_analysis/vol4/bpow/references/GCF_000001405.26_GRCh38_genomic.fna inputs/S3/L003/S3_L003_R1.fastq.gz inputs/S3/L003/S3_L003_R2.fastq.gz
@PG ID:MarkDuplicates-4B2F2425  VN:2.9.2-SNAPSHOT   CL:picard.sam.markduplicates.MarkDuplicates INPUT=[mapped/S3.L003.bam] OUTPUT=markdup/S3.L003.bam METRICS_FILE=markdup/S3.L003.metrics.txt REMOVE_DUPLICATES=false TMP_DIR=[/scratch/dfiler] CREATE_INDEX=true    MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 TAG_DUPLICATE_SET_MEMBERS=false REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag ASSUME_SORTED=false DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json    PN:MarkDuplicates
K00269:198:H55LWBBXY:2:2222:28392:34618 83  NC_000001.11    1180725 60  29S121M =   1180727 -119    ACTCTTTCCCTACACGACGCTCTTCCGATCTCCAGGAGAGCAGCTGCTGTACCAGCTTCCCAACAACAAGCTCCTCACCACCAAGATCGGGCTGCTCAGCACCCTTCGGGGACGGGCACGGGCCATGAGCAAGGCCAGCAAGGTGCCGGG  -JJJJJFJJJJJJJJJFJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA  MD:Z:121    NM:i:0  AS:i:121    XS:i:0  RG:Z:S3_L002    PG:Z:MarkDuplicates
K00269:198:H55LWBBXY:2:2222:28392:34618 163 NC_000001.11    1180727 60  119M31S =   1180725 119 CCAGGAGAGCAGCTGCTGTACCAGCTTCCCAACAACAAGCTCCTCACCACCAAGATCGGGCTGCTCAGCACCCTTCGGGGACGGGCACGGGCCATGAGCAAGGCCAGCAAGGTGCCGGGAGATCGGAAGAGCACACGTCTGAACTCCAGT  AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJ  MD:Z:119    NM:i:0  AS:i:119    XS:i:19 RG:Z:S3_L002    PG:Z:MarkDuplicates

## github.vcf ##
##fileformat=VCFv4.2
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  SMP1    SMP2
NC_000001.11    1180851 NC_000001.11:1180851_T/C    T   C   20450.4 .   AB=0.471229;ABP=12.2568;AC=1;AF=0.25;AN=4;AO=831;CIGAR=9M1X;DP=133;DPB=8557.7;DPRA=1.09321;EPP=98.3382;EPPR=1443.91;GTI=0;LEN=1;MEANALT=2.60976;MQM=60;MQMR=59.9988;NS=126;NUMALT=3;ODDS=13.1657;PAIRED=0.998797;PAIREDR=0.997442;PAO=430.75;PQA=14556.7;PQR=14751.7;PRO=442.083;QA=29918;QR=245335;RO=6646;RPL=612;RPP=406.598;RPPR=3656.66;RPR=219;RUN=1;SAF=117;SAP=934.337;SAR=714;SRF=742;SRP=8709.24;SRR=5904;TYPE=snp;technology.Illumina=1  GT:AD:AO:DP:GQ:PL:QA:QR:RO  0/0:76,0:0:78:99:0,436,2669:0:2983:76   0/1:26,27:27:55:99:594,0,542:1065:1007:26

$ variant github.sam -l github.vcf 
@HD VN:1.5  SO:coordinate
@SQ SN:NC_000001.11 LN:248956422
@RG ID:S3_L002  SM:S3   PL:Illumina CN:UNC
@RG ID:S3_L003  SM:S3   PL:Illumina CN:UNC
@PG ID:bwa  PN:bwa  VN:0.7.15-r1140 CL:bwa mem -t 8 -R @RG\tID:S3_L002\tSM:S3\tPL:Illumina\tCN:UNC /projects/sequence_analysis/vol4/bpow/references/GCF_000001405.26_GRCh38_genomic.fna inputs/S3/L002/S3_L002_R2.fastq.gz inputs/S3/L002/S3_L002_R1.fastq.gz
@PG ID:MarkDuplicates   VN:2.9.2-SNAPSHOT   CL:picard.sam.markduplicates.MarkDuplicates INPUT=[mapped/S3.L002.bam] OUTPUT=markdup/S3.L002.bam METRICS_FILE=markdup/S3.L002.metrics.txt REMOVE_DUPLICATES=false TMP_DIR=[/scratch/dfiler] CREATE_INDEX=true    MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 TAG_DUPLICATE_SET_MEMBERS=false REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag ASSUME_SORTED=false DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json    PN:MarkDuplicates
@PG ID:bwa-4A25094  PN:bwa  VN:0.7.15-r1140 CL:bwa mem -t 8 -R @RG\tID:S3_L003\tSM:S3\tPL:Illumina\tCN:UNC /projects/sequence_analysis/vol4/bpow/references/GCF_000001405.26_GRCh38_genomic.fna inputs/S3/L003/S3_L003_R1.fastq.gz inputs/S3/L003/S3_L003_R2.fastq.gz
@PG ID:MarkDuplicates-4B2F2425  VN:2.9.2-SNAPSHOT   CL:picard.sam.markduplicates.MarkDuplicates INPUT=[mapped/S3.L003.bam] OUTPUT=markdup/S3.L003.bam METRICS_FILE=markdup/S3.L003.metrics.txt REMOVE_DUPLICATES=false TMP_DIR=[/scratch/dfiler] CREATE_INDEX=true    MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 TAG_DUPLICATE_SET_MEMBERS=false REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag ASSUME_SORTED=false DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json    PN:MarkDuplicates
K00269:198:H55LWBBXY:2:2222:28392:34618 83  NC_000001.11    1180725 60  29S121M =   1180727 -119    ACTCTTTCCCTACACGACGCTCTTCCGATCTCCAGGAGAGCAGCTGCTGTACCAGCTTCCCAACAACAAGCTCCTCACCACCAAGATCGGGCTGCTCAGCACCCTTCGGGGACGGGCACGGGCCATGAGCAAGGCCAGCAAGGTGCCGGG  -JJJJJFJJJJJJJJJFJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA  MD:Z:121    NM:i:0  AS:i:121    XS:i:0  RG:Z:S3_L002    PG:Z:MarkDuplicates
K00269:198:H55LWBBXY:2:2222:28392:34618 163 NC_000001.11    1180727 60  119M31S =   1180725 119 CCAGGAGAGCAGCTGCTGTACCAGCTTCCCAACAACAAGCTCCTCACCACCAAGATCGGGCTGCTCAGCACCCTTCGGGGACGGGCACGGGCCATGAGCAAGGCCAGCAAGGTGCCGGGAGATCGGAAGAGCACACGTCTGAACTCCAGT  AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJ  MD:Z:119    NM:i:0  AS:i:119    XS:i:19 RG:Z:S3_L002    PG:Z:MarkDuplicates

For illustration, using ASCIIGenome:

$ ASCIIGenome github-out.sam github.vcf
github-out.sam@2; Reads: 2                                                                                                          
ctccaggagagcagctgctgtaccagcttcccaacaacaagctcctcaccaccaagatcgggctgctcagcacccttcggggacgggcacgggccatgagcaaggccagcaaggtgccggg           
  CCAGGAGAGCAGCTGCTGTACCAGCTTCCCAACAACAAGCTCCTCACCACCAAGATCGGGCTGCTCAGCACCCTTCGGGGACGGGCACGGGCCATGAGCAAGGCCAGCAAGGTGCCGGG           
github.vcf#3; N: 1                                                                                                            C     
SMP1                                                                                                                          .     
SMP2                                                                                                                          E     
1180725   1180735   1180745   1180755   1180765   1180775   1180785   1180795   1180805   1180815   1180825   1180835   1180845   11
0         .08       .15       .23       .31       .38       .46       .53       .61       .69       .76       .84       .92       .9
NC_000001.11:1180725-1180856; 132 bp; 1.0 bp/char                /\

You can see the mapped portion of the reads do not cover the locus of interest.