Will_Gilks ✭✭

About

Username
Will_Gilks
Location
University of Sussex, UK
Joined
Visits
443
Last Active
Roles
Member
Points
78
Badges
7
Location
University of Sussex, UK
Twitter
@https://twitter.com/wpgilks
Full Name
William Gilks

Comments

  • Hi @sysuls Is that ggplot you're using ? To obtain a better view of your data, firstly, format the GATK variants info table as 'long' using dplyr with the tidyr-gather function : long.df = df %>% gather (metric, value, -c(CHROM, POS, ID)) Then…
  • I'd just like to add that contrary to previous indications, no tears were shed as a result of Genomestrip (or even GATK); rather these implications were made in jest.
  • Working solution below. Collapsing coordinates of overlapping events is possible using bedtools but the regenotyping proved tricky. ## DESCRIPTION## Bash shell script for processing of vcf genotype files generated by Genomestrip/2.0## Specifically …
  • (Quote) Yes, I've attached a png and pdf of a heterozygous 4Kb duplication, identified and genotyped using the GS CNV pipeline (not the GS SV pipeline). Beware of single events that GS splits into multiple.
  • Hi, I was wondering what all the metrics stand for in the vcf files outputted by the Genomestrip SV and CNV pipelines. Apologies if this is in the documentation and I've missed it. I know some are fairly self-explanatory, but other such as cohFN, m…
  • You should use the same reference genome for both read-mapping and SNP-calling. Alternatively you can remove unmapped scaffolds from the current SNP-calling reference genome but I don't recommend it.
  • Thanks Sheila. This code seems to do what I need. ## Combine vcfs GenomeAnalysisTK -R ${refseq} \ -T CombineVariants \ -genotypeMergeOptions UNSORTED \ …
  • @sumedhagarg I use this code for making the various reference genome helper-files. Ideally the helper-files would only be made once by the same group that assembled the genome. This would prevent errors caused by people using different methods, and…
  • Hi @paumarc To calculate number of reads per chromosome - which will tell you indirectly the sex of your subjects - try Samtools idxstats http://samtools.sourceforge.net/ with a loop to run through the bam files. for i in *.bam; dosamtools idxstat…
  • Hi @bhandsaker I added the advised parity correction disable but the result wasn't noticeabley different. Code as ## CNV pipeline ____________________________________________ mkdir -p cnv_PAR_pipeline || exit 1 mkdir -p cnv_PAR_pipel…
  • Hi @bhandsaker GenerateHaploidCNV genotypes worked. I've attached screenshots from the 'Integrated Genomics Viewer' of one whole chromosome arm of the resulting vcf. As before, homozygous variants are expected to be very rare due to our breeding d…
  • Hi @bhandsaker Many thanks, and happy travels.
  • @bhandsaker I'm assuming you meant export SV_DIR=/cm/shared/apps/svtoolkit/2.0.1602/. Anyway, returns the "can't find R script" under any variation I try. Should I be able to locate the script in /cm/shared/apps/svtoolkit/2.0.1602/ ? It'…
  • Hi @bhandsaker Using -intervalsList for the five major scaffolds/chromosomes, the CNV calling pipeline completed in about 3 days. I've attached the variants-per-sample pdf report from stage 5 evaluation. There seems to be a few samples with elevate…
  • Hi @bhandsaker I chose to plot a selection of interesting looking events, mostly duplications, so it's not an unbiased selection. I haven't purposely disabled parity correction. To be honest, I haven't heard of it. Hardy-Weinberg Equilibrium is ex…
  • @bhandsaker Thanks. It worked with -intervalList test_intervals.list in which the .list file contained only "chr2L". The hanging at stage 12 is resolved. The only indication that the job was running at this stage was increase tmpdir/ fil…
  • Thanks, good points. Using GATK, I concatenated vcfs from each chromosome, then filtered by AN and AC. The END attributes are retained for normal variant records, just not this example. In the virgin vcf, the two records which later become merged ha…
  • @RLCollins I haven't tried any of the above yet. Is it compatible with non-human organisms? I suspect a bit more tweaking is required.
  • @Sheila I have a similar problem. Multiple of records with no MQRankSum or ReadPosRankSum value are passing filters. I see how records with missing values can be forced to fail filters by following the link previously provided (using --missingValu…
  • @adr1an I've found the name label for the mitochondrial genome can between assemblies. Eg chrM vs mitochrondrail_genome. Also reference meta-data name labels can vary within an assembly version. Eg between fasta and chain file types. Also, it's quit…
  • Somehow this post got marked as answered when it isn't ... :(
  • I came to a similar point but am now held up due to technical isssues, so best of luck to you. The response is here. http://gatkforums.broadinstitute.org/gatk/discussion/7208/genotyping-cnvs#latest I think you basically need to run GenerateHaploid…
  • As far as I can make out it means that at that position, there is neither a C or a CGTT - so basically there's a deletion of the C allele. It wouldn't surprise me if the vast majority of these kind of compound-null loci are generated in error - and …
  • A problem with the bam suggested at the end of the error message perhaps? ie. genstrip_chr12:1-133275309:100-100000_pairsbyname1_4806193241021307727.bam .. so run with about half the bam files, or check their format and quality. Or maybe there's a …
  • @cwhelan Thanks Chris, very informative. I guess the problem with trying to use -intervalList was due to not having a reference bundle, or only having custom masking and parameters files for the whole genome rather than for each chromosome .. but th…
  • Ta, could've sworn I put this in genomestrip.
  • To get a genotypes table with 0, 1 and 2 representing genotypes, I used Plink to generate a binary dosage format. https://www.cog-genomics.org/plink2/data ## Biallelic only variants ___ (removes null alternates hated by plink)GenomeAnalysisTK -R ${…
  • @bhandsaker My sample is unusual in that all individuals have a mother from the reference assembly strain, and different fathers from an wild, outbred population. DEL2224 is genotyped as heterozygous in the maternal (reference) line. It seems that t…
  • @bhandsaker Thanks for the tips. Actually, I was anticipating three plots, with the two others showing results from different methods, and read depth in individuals across the event as shown in http://www.broadinstitute.org/software/genomestrip/site…
  • @bhandsaker Partition file seems to me in place. I'm nearly there but aren't there meant to be two other graphs for this output? I've tried --pretty TRUE, -pretty true etc. Is the final vcf required at all here, as I customised the variant ids in it.
  • Hi @bhandsaker I've attached pdf and png of various deletion metrics across the genome. They've been filtered as: GSCOHERENCE > -100, GSDEPTHPVALUE < 1e-06. Not shown is a 300Kb event on chr3R. There's still a lot of events detected around t…
  • Behold, I have found a solution not two hours later, with the addition of --unsafe LENIENT_VCF_PROCESSING \
  • this post has been answered
  • This post has been answered ...
  • Whatever my problem was this code now werks: #!/bin/sh#$ -N z_gstrip_deletions#$ -pe openmp 30#$ -S /bin/sh#$ -cwd#$ -j y#$ -q bioinf.q## Genomic structural variation identification and genotyping.## Drosophila melanogaster n=221## WGS-NGS data, al…
  • Hi @Sheila @Geraldine_VdAuwera I have two questions: * The recommended hard-filters between SNPs and indels differ conspicuously. This is most prominent with the tenfold relaxation of Fisher Strand bias threshold. I suspect this is because of the …
  • @Michele I had that error message recently because I'd missed a new-line backslash. GATK tolerates 0/* but other programs don't.
  • @niranjanks 1. Get start and stop coordinates for all indels. 2. Expand coordinates using e.g. R and export as a .list 3. gatk filterVariants adding list as a mask I think, or exclude regions. 4. Gatk selectVariants
  • @bhandsaker Looking at the most common deletion passing filters, using IGV, it's pretty convincing, with greater distance between pairs (see attachment). I looked at a 10K deletion that failed filters, and is only present in one individual, and tha…
  • @bhandsaker ... maybe the attached pdf less blurry than the previous png ...?
  • @bhandsaker Thanks, all that's encouraging. The phasing of CNVs and SNPs is exciting. I'm looking through Beagle now to try and understand what's required. I'm also trying to submit the CNV and deletion records from Genomestrip, to NCBI dbVar but t…
  • @daveuu I solved some Jexl expression filter issues by removing the whitespace from the expressions themselves. Maybe this will help you. I think it's sensitive to linux version or something. Will.
  • ... solved by doing only the major chromosomes (not small unmapped scaffolds) in separate discovery runs.
  • Whatever the exact cause of this error, I've circumvented it by running Deletion Discovery on the major chromosomes separately. Assembly dm6 has a lot of small unmapped scaffolds which I guess are causing problems, perhaps if they are smaller than t…
  • Thanks, works with: GenomeAnalysisTK -R local_reference/dm6.fa \ -T VariantEval \ -eval ${myvcf} \ --evalModule CountVariants \ --stratificationModule Sample \ -noEV \ --out ${myvcf}.vareval.txt \ Output shown in image.
  • Using ... GenomeAnalysisTK -R local_reference/dm6.fa \ -T VariantEval \ -eval ${myvcf} \ --evalModule CountVariants \ --stratificationModule Sample \ --out ${myvcf}.vareval.txt \ I'm getting ... ERROR MESSAGE: Invalid command line: Argument ST and…
  • @laoyeye Did you ever find a solution to this problem? I'm passing -runDir. I've got combined coherence.dat, depth.dat etc files but the unfiltered.vcf is empty.
  • Thanks both, sounds great. I'll give it a try and post the result when I have a spare.
  • What ever my problem was there, I've fixed it, so the question can be marked as answered.
  • Thanks @Geraldine_VdAuwera If I understand you correctly, one of the peaks should be full of variants which are mostly heterozygous, and the other peak is full of variants that are mostly homozygous. I shall try and plot the quality distribution by …
  • @says_anova My problem was mostly down to typos. Check the input vcf (or vcf list) is complete, and code that generated it is correct. I was genotyping gvcfs but it seems you're doing something different.
  • Hallelujah. Made ploidy file as you advised. svmask required an indexed fasta but easy to do. Below is some code for generating the svmask on a linux SGE without the hashed headers in case anyone else needs it. module load sge module load genomest…
  • Thanks @bhandsaker Have corrected ploidy file and made rdmask. To make the svmask do you still recommend computing the mask for each chromosome separately, and then concatenate? (Each of the four major autosomes are about 30MB.) http://gatkforums.b…
  • The other thing I was wondering about was Samtools. I read somewhere that GS requires it but I don't see it specified in any command line. I'm currently using v1.2.
  • Deleted other files in the reference assembly directory but still similar error messages. Starts with an R error in preprocessing. Not sure if this is critical though: INFO 17:00:47,638 QJobsReporter - Plotting JobLogging GATKReport to file /lustre…
  • Thanks for getting back to me. I've attached the suggested metadata. There are other files stored in the same directory as dm6.fa. I'll delete these out as they're not for GS, and let you know how things go.
  • Did the negative unobsOut problem ever get solved ? I have a similar problem myself.
  • @bhandsaker Thanks for getting back to me. * Got it but have other stuff to fix. (Deletion pipeline atm. 20-30X, all females) * Got it. (Reads are long-ish - will exclude alt allele stage). * Got it. (Preprocessing on as much of genome as possi…
  • To answer my own question: To make a vcf of variants in repeats only: GenomeAnalysisTK -R local_reference/dm6.fa \ -T SelectVariants \ -V somones.vcf \ -L ../masking_intervals/dmel6_repMask.gatk.intervals \ -o repeats.vcf To mark a…
  • @Sheila Ok, thanks for the link, I must've missed before. Example code below. GenomeAnalysisTK -R ${ref_genome} \ -T VariantFiltration \ -V raw.vcf \ -filter "QUAL<1000.0" -filterName "lowQual&q…
  • Your paths are rather long. Could easily be a mistake in them. Is that really all of the information in the error message? Do the numbers refer to line numbers in the vcf, which might indicate the origin of the error.
  • * dbSNP have informed me that they don't know what <*> stands for. dbSNP currently doesn't accept any variants with non-ATGC in either ref. or alt. alleles .... which makes it clear. I guess just delete those variants then. * The simple way …
  • Hi @Sheila Joy but not nirvana. DONE: Excluding by event length less than 50bp: GenomeAnalysisTK -R ${refgenome} -T SelectVariants -V ${myinvcf} -o nolong.vcf --maxIndelSize 50 N.B. Official documentation incorrect with '-maxIndelSize'. Suggest co…
  • Thanks @Sheila clearly a matter of me looking in the correct places. GATK guides seem to be becoming rather scattered. How about my question 2.. How to set individual genotype calls which have less than 10X coverage to ./. ?
  • Naming variants in a vcf file My bash solution below. Would be interested if anyone has a better way. Naming variants by position is not necessarily a great idea as they can change by build version, so alternatively just make a list of unique rando…
  • Hi @Sheila Well, really I have two problems in one. My first problem is lifting over a vcf to a more recent reference. I've solved this, using a combination of GATK (LiftOverVariants and FilterLiftedVariants), vcfsorter.pl, bash and reading into py…
  • Hi, I'm trying to convert the D.melanogaster genotype data from the DGRP collection from dm3 to dm6. I've run the perl script to sort the vcf as the reference but I'm 'getting chain file not compatible'. Chain file is dm3ToDm6.over.chain obtained f…
  • The problem in this case was space. Bad: --filterExpression "QUAL < 1000.0" --filterName "LowQual" Good: --filterExpression "QUAL<1000.0" --filterName "LowQual"
  • Hi @Sheila Thanks for your endeavours. I have a problem using gatk filter expressions on the SGE cluster but not my desktop. (See link) I've pretty much excluded every other cause. I'll go back to it in a bit. http://gatkforums.broadinstitute.org/di…
  • @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 genera…
  • 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…
  • 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.
  • 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 val…
  • I'm using gatk/3.4-46 on both cluster and desktop. The above command and error message are from the cluster but it generally works fine on my desktop. I removed the speech-marks (") from the -G_filter "isNoCall > 0.0" which change…
  • @annat @baesc @drramki I've got an incomplete VCF header problem as well (using gatk v3.4-46): Program Args: --analysis_type CombineGVCFs --reference_sequence ../../reference_sequences/dmel/v6.0/dm6.fa --variant rg_GVCF.list Returns: Key END foun…
  • Thanks @thibault . I think that was our main problem. Successful input shown below: module load sge module load genomestrip/2.0 java \ -cp /cm/shared/apps/svtoolkit/2.0.1602/lib/SVToolkit.jar:/cm/shared/apps/svtoolkit/2.0.1602/lib/gatk/GenomeAnaly…
  • Sorry, I don't think I should've tagged this one as "Answered". Any chance you could change this?
  • Hi @Sheila Thanks for your response. In fact I haven't performed BQSR yet. It seems like a good thing to do. I trimmed reads (prior to mapping with BWA-mem and then Stampy) with: fastq-mcf \ ~/adaptor_sequences/adaptors_truseq_130214.fast…
  • The two big spikes in the error-rate-by-cycle plot, lower-left (E): Variants in reads generated during these cycles are more likely to be spurious, right? So the genotypes should be filtered out - or the sample discarded. There's a few other samples…