malariagen / pipelines Goto Github PK
View Code? Open in Web Editor NEWPipelines for processing malaria parasite and mosquito genome sequence data.
License: MIT License
Pipelines for processing malaria parasite and mosquito genome sequence data.
License: MIT License
We need to reduce the peak size of the vector alignment pipeline so that we can run more pipelines in parallel without running out of dis. Peak size refers to the total size of all files generated by the pipeline, including intermediate files. There is built-in way for cromwell to delete intermediate files as we go along. So by default, we keep all intermediate files until a workflow is complete.
One way to reduce peak size is
Verify that the output crams are lossless with respect to information from the bams. Verify that genotyping pipeline is not affected by use of crams in alignment pipeline.
The pertinent files are anything related to alignment. The possible workflow WDL entry points.
Follow the entrypoint WDLs to find all the dependent WDLs that may need modifying.
See https://github.com/malariagen/pipelines/blob/master/docs/specs/short-read-alignment-vector.md for alignment pipeline specs.
The mandatory fields in the manifest are sample_id, irods_path. However, their order should not matter. Also, if a user decides to add extra fields for their own use, it should not matter. From comments from pull request: #67 (comment)
We need to solidify the final manifest fields required for vector alignment and genotyping. Right now, we are using fields found in https://github.com/malariagen/pipelines/blob/master/pipelines/short-read-alignment-and-genotyping-vector/farm5/input_files/validation/validation_lanelets_farm5_local.tsv
We don't need all these fields for the production manifest passed in to the batch-irods import-alignment-genotyping vector pipelines. At the moment, we probably only need:
However, run_ena is used as an input in the per-lanelet-alignment vector pipelines.
We will need to solidify the final inputs in the manifest that we need to keep track of samples and their progress in the production system. For example, we don't need the study in the input manifest to run the pipelines. But will we need the study stored as a workflow label/metadata so that we can query for it later in the future?
Share example mosquito crams released from Core Sequencing with @gbggrant and @kbergin.
If possible, kill birds with one stone and try to find the smallest valid crams possible to solve #25 . Querying and sorting IRODS by file size may be a bit fiddly. Instead, let's choose crams with the fewest reads according to sample metadata from past releases.
@gbggrant requested ways to make bwa short read alignment deterministic for testing.
Trawling through the internet indicates that we can make bwa deterministic if we follow all of these steps:
Always use the same fastq with the reads in the same order
Always use the same machine and bwa settings. Alternatively, if we need to change the machine and threads given to bwa, we can explicitly set the hidden bwa mem option -K to set the chunk size. If -K is not set, then the chunksize is determined by the number threads. This makes a difference to the end bam because the insert size is determined by the chunk.
See past comments from Heng Li:
lh3/bwa#192
lh3/bwa#272
lh3/bwa#121
Also see this centre for common disease genomics repo issue:
CCDG/Pipeline-Standardization#2
For the mosquito short read alignment pipeline, I believe there are two main options regarding the input data source for the sequence reads:
There are some practical trade-off in these choices:
Copy script from https://github.com/malariagen/legacy_pipelines/blob/master/prod-tools/scripts/vcf_to_zarr.py to this repository.
@gbggrant All the legacy scripts used the by vrpipe (the legacy workflow management system) should be found in https://github.com/malariagen/legacy_pipelines in case you want to refer to them.
The legacy vcf_to_zarr.py script needs edits for the new vector pipelines, which we will review in a PR:
Part of #2 but breaking this out as a separate issue for convenience. Identify a suitable dataset to use for testing and validation purposes for the short read alignment pipeline.
This is a placeholder for discussion how to version workflows.
Given a batch sample manifest file with multiple samples, each row specifying lanelet irods path and sample ID, download all sample lanelets in parallel from IRODS, and run each sample through alignment and genotyping pipeline.
The initial version of the pipeline generated an output bam that failed ValidateSamFile due to a difference between the implementations of the calculation of the value of the NM tag in samtool's calmd tool and that in picard's ValidateSamFile tool. See: samtools/samtools#717 (comment)
After discussion with @alimanfoo and others, and with the understanding that the NM tag is not used elsewhere (particularly in gatk 3's UnifiedGenotyper), the decision was made to remove the usage of calmd and replace it with Picard's SetNmMdAndUqTags tool.
Suggestions from @jon.keatley.aspis on slack Mar 26
Jon Keatley 16:17
ssh farm5-login
ls /lustre
check for .irods in ~/
busers & bugroup
Any new user should have access to:
Add dockerfiles for all the tools needed by the vector short read alignment pipeline version 1.1.1.
Part of #2, breaking out as separate issue for convenience.
Add some kind of CI validation of WDL files to check syntax etc.
Legacy pipelines overwrite the read group ID, library, study, and sample name from the released bams/crams from Sanger Core Sequencing. Determine which fields we want to overwrite, going forward, if any.
We anticipate that this repository will store all code artefacts necessary for defining data pipelines. E.g., WDL files, docker files, test datasets, pipeline documentation, etc. What file and folder structure should we use for organising these? Any other conventions we need to agree?
Currently the pipeline wdl files for gcp and farm5 are extremely similar. It would be great to make these one file and pass an argument to the build_pipeline_release script which would specify the platform to generate the wdl for.
A concern about this is that in making the pipeline wdl platform non-specific (for instance by adding a placeholder field for the platform-to-be-specified) this might make the wdl non-valid for our wdl checking tool (miniwdl).
Dear MalariaGen Team
Apologies, if I am raising this question in wrong repository. I tried to download P. vivax variant data (PV4 Dataset) from https://www.malariagen.net/resource/30. I am actually interested in using this data for Base Quality Score Recalibration (BQSR) step of GATK germline variant calling pipeline. I tried to merge the Per-chromosome variant files and then I filtered the merged variant file for high confidence variants (SNPs and Indels) as recommended in the README file. The bcftools stats
shows correct number to total variants and samples. However, the Number of SNPs reported in bcftools stats report are different than what has been reported in README. Can you help me understand why is it so. I am also sharing the code used below
###############Bcftools ####################
SN [2]id [3]key [4]value
SN 0 number of samples: 1895
SN 0 number of records: 1303984
SN 0 number of no-ALTs: 0
SN 0 number of SNPs: 1117804
SN 0 number of MNPs: 0
SN 0 number of indels: 254148
SN 0 number of others: 0
SN 0 number of multiallelic sites: 331817
SN 0 number of multiallelic SNP sites: 33748
##############################README
The VCF files contains details of 4,571,056 discovered variant genome positions.
These variants were discovered amongst all samples from the release.
3,083,454 of these variant positions are SNPs, with the remainder being either
short insertion/deletions (indels), or a combination of SNPs and indels. It is
important to note that many of these variants are considered low quality. Only
the variants for which the FILTER column is set to PASS should be considered of
reasonable quality. There are 1,303,984 such PASS variants of which 945,649
are SNPs and 358,335 indels.
wget ftp://ngs.sanger.ac.uk/production/malaria/Resource/30/Pv4_vcf/*
ls -1 *.vcf.gz > vcf.list
#merging vcfs
picard GatherVcfs I=vcf.list O=merged.vcf.gz
gatk IndexFeatureFile -I merged.vcf.gz
# using only high confidence variants for BQSR
gatk SelectVariants -R ../results/00_indexes/bwaidx/PlasmoDB-50_PvivaxP01_Genome.fasta -V merged.vcf.gz --exclude-filtered true -O filter.vcf.gz
## Complete output is attached below
bcftools stats filter.vcf.gz > filter.bcftoolstats.txt
Assembly used: PlasmoDB-50_PvivaxP01_Genome.fasta from PlasmoDB
gatk4-4.2.6.1-1
bcftools: Version: 1.10.2 (using htslib 1.10.2-3)
filter.bcftoolstats.txt
This issue is intended to cover implementation of a pipeline for producing analysis-ready alignments from paired-end short read Illumina sequencing of Anopheles moquitoes.
Some tasks (not necessarily exhaustive):
from @gbggrant
I've run the 167 Burkina Faso samples through the phasing pipeline. I haven't yet integrated your cohort_vcf_to_zarr script, but the VCFs for 2R and 3R are at:
gs://malariagen/Phasing/outputs/Ag1000Phase2_BurkinaFaso/Ag1000Phase2_BurkinaFaso_2R_phased.vcf.gz and
gs://malariagen/Phasing/outputs/Ag1000Phase2_BurkinaFaso/Ag1000Phase2_BurkinaFaso_3R_phased.vcf.gz - full outputs for the pipeline can be found at gs://malariagen/Phasing/outputs/Ag1000Phase2_BurkinaFaso/ - let me know if you can't read these (I can download them to lustre).
@hardingnj to run H12 scan, and check we see signals of natural selection we expect.
Issue to track production run of the mosquito phasing pipeline to build haplotypes for the Ag 3.0 (a.k.a. Ag1000G phase 3) data release.
Check that genotypes from new WDL pipelines are >99% similar to genotypes from legacy pipelines for the same samples.
Legacy pipeline zarrs can be found on Sanger lustre:
/lustre/scratch118/malaria/team112/pipelines/setups/vo_agam/symlink/vo_agam_vcf2zarr/convert_vcf_to_zarr_serial_contig
Scripts that check genotype concordance for parasite release zarrs can be found: https://gitlab.com/malariagen/gsp/pf7.0/-/blob/master/work/59_genotype_concordance_script/genotype_concordance.py
Note parasite release zarrs use a different schema from mosquito release zarrs. At the moment, the script will crash when it runs on mosquito zarrs.
@gbggrant suggested we provide a minimal input so that we can run the vector short read alignment pipeline in under an hour. At present, the small bams the malaria programme has provided will process in at least 8 hours on the Broad google cloud infrastructure.
The minimal input need not be scientifically valid. They only need to test that the pipelines are connected properly.
For legacy pipelines, we tested plumbing by using small sample test bams. I would share those, but the legacy pipelines were originally developed for human sequencing, and so the test bams were actually downsampled human bams. For these pipelines, let's try downsampling the mosquito crams.
I suspect that randomly downsampled crams can easily be used to test short read alignment and genotyping. However, we may need to ensure there are sufficient genotypes after genotyping to test phasing without crashing the phasing software.
Fix so that it only has a single .bam suffix
When running against an initial 28 samples from the actual AG3.2 dataset (https://github.com/malariagen/vector-ops/blob/master/tracking/1177-VO-ML-LEHMANN-VMF00015/wgs_lanelets.tsv), we noticed that the disk usage was > 300GB / sample.
Part of this is because files from one task are copied over as inputs for the next task, doubling the disk usage.
Change cromwell server config such that localization = soft-link. Check that the input files are symlinked from the output of the previous task. Check that this happens for all tasks, including the first task. In the past, we've noticed that the WDL will copy over the initial WDL input files, regardless of the localization config.
Add the changes made to the Picard Dockerfile by @gbggrant to resolve the location issue when using a singularity container.
Picard was stored by default in /root
but because singularity executes with the user's credentials it is impossible to access unless executed as root. The change proposed by @gbggrant is to move the jar from /root
to /bin
Umbrella issue to communicate re development of the mosquito phasing pipeline. Feel free to raise new issues to deal with specific topics and link back here.
The goal is to establish a baseline pipeline for comparison by porting https://gitlab.com/malariagen-aspis/aspis-pipeline/-/tree/master/Stage1-Pipeline with minimal functional changes.
Note that "Stage 1" encompasses BCL -> VCF and is broken into 3 steps, roughly: 1) BCL -> unaligned CRAMs, 2) -> aligned CRAMs, 3) -> pileup-based SNP calling using BCFtools. Methods improvements will focus on step 3, so this initial pipeline will only consist of this step.
It is likely that the result will be a single-sample WDL for step 3. We can file a subsequent issue to rework steps 1 and 2 as necessary.
What datasets are available for calibrating, evaluating and/or comparing different methods for variant calling and/or filtering, in Plasmodium or Anopheles?
Issue for discussion of work to compare outputs of the short read alignment pipeline with expected outputs from previous pipeline implementation.
On Sanger farm5, use LSF group pathdev instead of malaria-dk for production
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.