Attention:
The frontline support team will be unavailable to answer questions until May27th 2019. We will be back soon after. Thank you for your patience and we apologize for any inconvenience!

remove any variant with a ./.

Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

Does anyone know how to remove a whole variant record in which the call for any of the samples is ./. ?

It should be simple but have been trying filter-select and select+jexl for a while now but to no avail.

Best Answers

Answers

  • pdexheimerpdexheimer Member ✭✭✭✭

    Why don't you select on AN? If you only select the variants where AN=(ploidy*nSamples), it will only grab the records where all of the samples are called

  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    Thanks @pdexheimer Filtering by AN is a good hard-coded solution.

    Unfortunately my code:
    GenomeAnalysisTK -R local_reference/dm6.fa -T SelectVariants -V gd_Samps_raw_SNPs.vcf -select 'AN = 438' -o aajexeltest.vcf

    Is returning:
    Invalid argument value '=' at position 8.
    Invalid argument value '438' at position 9.

    ..whether or not I use one 'equals' sign or two. I'm not sure if I'm missing something.

  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    Also the other related query I have is how to remove any variant records than have a 1/1 call for any of the samples.... It would have to a different command to altogether.

  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    Filtering by AN now works well (on desktop) with:

    java -jar -Xmx2g ~/SOFTWARE/GenomeAnalysisTK-3.4-46/GenomeAnalysisTK.jar -R ../reference_sequences/dmel/v6.0/dm6.fa -T VariantFiltration -V gd_Samps_raw_SNPs.vcf --filterExpression "AN < 438" --filterName "lowCallRate" -o aa_marked_test.vcf

    I'll try on the cluster in a minute. Suspect you are correct about command syntax error. Not trying to call from another program, just running as a job in SGE.

    Thanks ever so much for the help.

    I don't suppose you have an answer for my second problem ? i.e. filtering variants that have a 1/1 call in an sample... I've done this in bash but I think it creates problems with the .idx

  • pdexheimerpdexheimer Member ✭✭✭✭

    It will cause problems for the .idx, yes - but you can just delete that and allow GATK to regenerate it next time you need it.

    The only way I know to do this in GATK is to iterate over every single sample, which is not something you want to do with 219 samples (--filterExpression 'vc.getGenotype("Sample 1").isHomAlt() | vc.getGenotype("Sample 2").isHomAlt() | ... - that's from memory, so the syntax may be a bit off, but you get the idea). With that many samples, you may be better off just using grep

  • KurtKurt Member ✭✭✭

    I found that when running --filterExpression from a shell script (like submitting to SGE). I have to replace double quotes with single quotes in my environment to get it to work.

    $JAVA_1_7/java -jar $GATK_DIR/GenomeAnalysisTK.jar \ -T VariantFiltration \ -R $REF_GENOME \ --variant $CORE_PATH/$PROJECT/TEMP/$SM_TAG"_QC_RAW_OnBait_SNV_ANNOTATED.vcf" \ --filterExpression 'QD < 2.0' \ --filterName 'QD' \ --filterExpression 'MQ < 30.0' \ --filterName 'MQ' \ --filterExpression 'FS > 40.0' \ --filterName 'FS' \ --filterExpression 'MQRankSum < -12.5' \ --filterName 'MQRankSum' \ --filterExpression 'ReadPosRankSum < -8.0' \ --filterName 'ReadPosRankSum' \ --filterExpression 'DP < 8.0' \ --filterName 'DP' \ --logging_level ERROR \ -o $CORE_PATH/$PROJECT/TEMP/$SM_TAG"_QC_FILTERED_OnBait_SNV.vcf"

  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    @kurt @pdexheimer Both ['] and ["] work fine for filter expression, from my mac, but neither are accepted on the linux SGE cluster. I've run other gatk commands in the same script prior to this error, and gatk is the only program used to generate the vcf files. I've tried moving the filtration command up the work-flow, in case the input vcf was malformed. I've re-typed GenomeAnalysisTK, checked the shell script headers and module environment, and made a shorter version of the whole script. I've simplified in the input file names.

    #!/bin/sh
    $ -N dummy_gatk
    $ -pe openmp 1
    $ -S /bin/sh
    $ -cwd
    $ -j y
    $ -q bioinf.q
    . /etc/profile.d/modules.sh
    module load jre/1.7.0_25
    module load gatk/3.4-0
    GenomeAnalysisTK -R ./local_reference/dm6.fa \
    -T SelectVariants -V comb_SNP_indel.vcf \
    -restrictAllelesTo BIALLELIC \
    -XL dmel6_repMask.gatk.intervals \
    -o bi_gd_Samps.vcf

    GenomeAnalysisTK -R ./local_reference/dm6.fa \ ##### problem line
    -T VariantFiltration -V bi_gd_Samps.vcf \
    --filterExpression 'AN < 438' --filterName 'lowCallRate' \
    -o bi_gd_Samps_lowCall_marked.vcf

    In the log is that there are no INFO-time-HelpFormatter rows before the error lines start which seems like a gatk loading problem (normal for previous and subsequent commands). I've asked sysadmin if there's anything odd regarding expressions on SGE.

    A USER ERROR has occurred (version 3.4-46-gbc02625):
    ##### ERROR MESSAGE: Invalid argument value '<' at position 8.
    ##### ERROR Invalid argument value '438' at position 9.

  • tommycarstensentommycarstensen United KingdomMember ✭✭✭
    edited August 2015

    @Will_Gilks There is a tool called bcftools, which can do this:

    bcftools view -g ^miss Will_Gilks.vcf.gz
    
  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    The problem in this case was space.

    Bad: --filterExpression "QUAL < 1000.0" --filterName "LowQual"

    Good: --filterExpression "QUAL<1000.0" --filterName "LowQual"

Sign In or Register to comment.