|
| 1 | +# Tutorial: Running a DivBase query on a public dataset |
| 2 | + |
| 3 | +This tutorial assumes that you have an account on DivBase and have a membership in a DivBase project with at least an EDIT role (i.e. can upload files and run queries). |
| 4 | + |
| 5 | +We will use a mouse (_Mus musculus_) data set availalbe on the European Nucleotide Archive: <https://www.ebi.ac.uk/ena/browser/view/ERZ022025>. It a 5.5 Gb VCF.gz file that contains 18 samples and 66,007,044 variants. |
| 6 | + |
| 7 | +## 1. Obtain the data and upload it to your DivBase project |
| 8 | + |
| 9 | +First we need download the files from EVA to our local computer. We can do that from the ENA homepage with a web browser, or from the terminal with a tool like `curl` or `wget`: |
| 10 | + |
| 11 | +```bash |
| 12 | +curl -o mgp.v3.snps.rsIDdbSNPv137.vcf.gz ftp://ftp.sra.ebi.ac.uk/vol1/analysis/ERZ022/ERZ022025/mgp.v3.snps.rsIDdbSNPv137.vcf.gz |
| 13 | +``` |
| 14 | + |
| 15 | +On the terminal, log in to DivBase: |
| 16 | + |
| 17 | +```bash |
| 18 | +divbase-cli auth login <USER_EMAIL> |
| 19 | +``` |
| 20 | + |
| 21 | +We then need to upload the file to the DivBase project. If the project you want to upload to is your default project, you can skip `--project <YOUR_DIVBASE_PROJECT_NAME>`, but otherwise you will need to specify that. |
| 22 | + |
| 23 | +```bash |
| 24 | +divbase-cli files upload mgp.v3.snps.rsIDdbSNPv137.vcf.gz --project <YOUR_DIVBASE_PROJECT_NAME> |
| 25 | +``` |
| 26 | + |
| 27 | +For the sake of demonstration, we can create a sidecar metadata file for the 18 samples in this dataset so that we can use that for the query. In fact, a file prepared for this demo can be downloaded from the DivBase repo with: |
| 28 | + |
| 29 | +```bash |
| 30 | +curl -o tutorial_mock_metadata_mgpv3snps.tsv https://raw.githubusercontent.com/ScilifelabDataCentre/divbase/refs/heads/docs-on-queries/tests/fixtures/tutorial_mock_metadata_mgpv3snps.tsv |
| 31 | +``` |
| 32 | + |
| 33 | +TODO update the TSV url when it has been merged to main. |
| 34 | + |
| 35 | +And then upload it to the DivBase project with: |
| 36 | + |
| 37 | +```bash |
| 38 | +divbase-cli files upload tutorial_mock_metadata_mgpv3snps.tsv --project <YOUR_DIVBASE_PROJECT_NAME> |
| 39 | +``` |
| 40 | + |
| 41 | +With this, we are set in terms of the files we need in the DivBase project to be able to run the query. |
| 42 | + |
| 43 | +### A comment on the mock sidecar metadata TSV |
| 44 | + |
| 45 | +If we look at the first few lines of this TSV, we see that the first column has the SAMPLE IDs from the VCF. There is a column for population grouping, and a column for the geographical area the sample comes from. All of this is made up for the tutorial, just to illustrate what a sidecar metadata file can look like. We will use the Area column later on to subset on samples. |
| 46 | + |
| 47 | +```bash |
| 48 | +head tests/fixtures/tutorial_mock_metadata_mgpv3snps.tsv |
| 49 | + |
| 50 | +#Sample_ID Population Area |
| 51 | +129P2 1 North |
| 52 | +129S1 2 East |
| 53 | +129S5 3 South |
| 54 | +AJ 4 West |
| 55 | +``` |
| 56 | + |
| 57 | +For comparison, we can inspect all the samples in the VCF with: |
| 58 | + |
| 59 | +```bash |
| 60 | +# On MacOS, use gzcat instead of zcat |
| 61 | +zcat mgp.v3.snps.rsIDdbSNPv137.vcf.gz | grep "#CHROM" |
| 62 | + |
| 63 | +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 129P2 129S1 129S5 AJ AKRJ BALBcJ C3HHeJ C57BL6NJ CASTEiJ CBAJ DBA2J FVBNJ LPJ NODShiLtJ NZOHlLtJ PWKPhJ SPRETEiJ WSBEiJ |
| 64 | +``` |
| 65 | + |
| 66 | +## 2. Submit a job to cache the VCF dimensions in the DivBase backend |
| 67 | + |
| 68 | +DivBase needs this to make calculations. It can be seen as a pre-processing step that needs to be run every time a new version of a VCF is uploaded to DivBase. |
| 69 | + |
| 70 | +```bash |
| 71 | +divbase-cli dimensions update --project <YOUR_DIVBASE_PROJECT_NAME> |
| 72 | +``` |
| 73 | + |
| 74 | +The task should now have been submitted. The terminal prints a DivBase Task ID with a message like this: |
| 75 | + |
| 76 | +``` |
| 77 | +# Example with Task ID 102 |
| 78 | +Job submitted successfully with task id: 102 |
| 79 | +``` |
| 80 | + |
| 81 | +Make note of the Task ID integer for for now; we will use it later in the tutorial for a variable named `<THE_JOB_ID_OF_THE_QUERY>`. |
| 82 | + |
| 83 | +We need to wait until this has finished before we can send the actual DivBase query. |
| 84 | + |
| 85 | +```bash |
| 86 | +divbase-cli task-history user |
| 87 | +``` |
| 88 | + |
| 89 | +## 3. Submit a query job |
| 90 | + |
| 91 | +```bash |
| 92 | +divbase-cli query bcftools-pipe --tsv-filter 'Area:North,East' --command 'view -s SAMPLES; view -r 1:15000000-25000000' --metadata-tsv-name tutorial_mock_metadata_mgpv3snps.tsv --project <YOUR_DIVBASE_PROJECT_NAME> |
| 93 | +``` |
| 94 | + |
| 95 | +note that we are specifically using the `tutorial_mock_metadata_mgpv3snps.tsv` sample metadata. If this is not specified, the query will default to `sample_metadata.tsv`. |
| 96 | + |
| 97 | +Depending on queue length, this job will take some time to run. Check that it has started, and then leave it running for half an hour. |
| 98 | + |
| 99 | +```bash |
| 100 | +divbase-cli task-history user |
| 101 | +``` |
| 102 | + |
| 103 | +## 4. Download the results file |
| 104 | + |
| 105 | +```bash |
| 106 | +divbase-cli file download result_of_job_<THE_JOB_ID_OF_THE_QUERY>.vcf.gz--project <YOUR_DIVBASE_PROJECT_NAME> |
| 107 | +``` |
| 108 | + |
| 109 | +We can now run some quick sanity-checks on the result file. |
| 110 | + |
| 111 | +```bash |
| 112 | +# On MacOS, use gzcat instead of zcat |
| 113 | +zcat result_of_job_<THE_JOB_ID_OF_THE_QUERY>.vcf.gz | grep -v "^#" |wc -l |
| 114 | + |
| 115 | +# Expected terminal output: |
| 116 | +297415 |
| 117 | +``` |
| 118 | + |
| 119 | +If you want, you can compare this to the original file: |
| 120 | + |
| 121 | +```bash |
| 122 | +# Note! This will take a little time since this file has many rows |
| 123 | +zzcat mgp.v3.snps.rsIDdbSNPv137.vcf.gz | grep -v "^#" |wc -l |
| 124 | + |
| 125 | +# Expected terminal output: |
| 126 | +66007044 |
| 127 | +``` |
| 128 | + |
| 129 | +As for the samples, the expected result is the following: |
| 130 | + |
| 131 | +```bash |
| 132 | +zcat mgp.v3.snps.rsIDdbSNPv137.vcf.gz | grep "#CHROM" |
| 133 | +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 129P2 129S1 129S5 AJ AKRJ BALBcJ C3HHeJ C57BL6NJ CASTEiJ CBAJ DBA2J FVBNJ LPJ NODShiLtJ NZOHlLtJ PWKPhJ SPRETEiJ WSBEiJ |
| 134 | + |
| 135 | +zcat result_of_job_<THE_JOB_ID_OF_THE_QUERY>.vcf.gz | grep "#CHROM" |
| 136 | +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 129P2 129S1 AKRJ BALBcJ CASTEiJ CBAJ LPJ NODShiLtJ SPRETEiJ WSBEiJ |
| 137 | +``` |
0 commit comments