# Question
Title: Minfi and estimateCellCounts Error: in p[trainingProbes, ]: subscript out of bounds
I would like to plot using estimateCellCounts as I do not know my cell counts for the blood samples I am using however, it seem that I cannot get estimateCellCounts to work with EPIC arrays.

Detailing Attempts:

**Attempt 1**:

* I first attempted it with EPIC arrays but it wouldn't take it.

```
estimateCellCounts(RGset,
                   compositeCellType = "Blood",
                   processMethod = "auto",
                   probeSelect = "auto",
                   cellTypes = c("Bcell",
                                 "CD4T",
                                 "CD8T",
                                 "Eos",
                                 "Gran",
                                 "Mono",
                                 "Neu",
                                 "NK"
                                ),
                   referencePlatform = "IlluminaHumanMethylationEPIC",
                   returnAll = TRUE,
                   meanPlot = TRUE,
                   verbose = TRUE)

```

And the associated error:

```
Loading required package: FlowSorted.Blood.EPIC
Loading required package: ExperimentHub
Loading required package: AnnotationHub
Loading required package: BiocFileCache
Loading required package: dbplyr

Attaching package: 'AnnotationHub'

The following object is masked from 'package:Biobase':
    cache

Warning message in data(list = referencePkg):
"data set 'FlowSorted.Blood.EPIC' not found"
Error in get(referencePkg): object 'FlowSorted.Blood.EPIC' not found
Traceback:
1. estimateCellCounts(RGset, compositeCellType = "Blood", processMethod = "auto",
 .     probeSelect = "auto", cellTypes = c("Bcell", "CD4T", "CD8T",
 .         "Eos", "Gran", "Mono", "Neu", "NK"), referencePlatform = "IlluminaHumanMethylationEPIC",
 .     returnAll = TRUE, meanPlot = TRUE, verbose = TRUE)
2. get(referencePkg)
3. get(referencePkg)

```

Adding the `library("FlowSorted.Blood.EPIC")` did not stop the error `"data set 'FlowSorted.Blood.EPIC' not found"`.

**Attempt 2**:

1. Then, I looked up the error I got first on Github and found that converting to 450k was recommended. So I used the `convertArray`.
2. Then it threw up an error saying `"FlowSorted.Blood.450k"` was not installed.

**Attempt 3**:

So I installed and loaded the library and re-ran to get the code below and the error below. However, this does not change the issue that it does not seem to run. The code I used and the error is below - any suggestions are appreciated:

```
library("FlowSorted.Blood.450k")

RGset_450k = convertArray(RGset,
                          outType = c("IlluminaHumanMethylation450k"),
                          verbose = TRUE
                         )

estimateCellCounts(RGset_450k,
                   compositeCellType = "Blood",
                   processMethod = "auto",
                   probeSelect = "auto",
                   cellTypes = c("Bcell",
                                 "CD4T",
                                 "CD8T",
                                 "Eos",
                                 "Gran",
                                 "Mono",
                                 "Neu",
                                 "NK"
                                ),
                   referencePlatform = "IlluminaHumanMethylation450k",
                   returnAll = TRUE,
                   meanPlot = TRUE,
                   verbose = TRUE)

```

**Error**:

```
Error in p[trainingProbes, ]: subscript out of bounds
Traceback:

1. estimateCellCounts(RGset_450k, compositeCellType = "Blood", processMethod = "auto",
 .     probeSelect = "auto", cellTypes = c("Bcell", "CD4T", "CD8T",
 .         "Eos", "Gran", "Mono", "Neu", "NK"), referencePlatform = "IlluminaHumanMethylation450k",
 .     returnAll = TRUE, meanPlot = TRUE, verbose = TRUE)

```

How would I make it get around this error as it does not output any results? Any suggestions/solutions are appreciated.

# Answer
I am the developer of FlowSorted.Blood.EPIC, I will try to provide you with some guidance. First, the FlowSorted.Blood.EPIC is located on ExperimentHub, which means that you need to ask the program to download the database first before running estimateCellCounts (minfi) or estimateCellCounts2 (FlowSorted.Blood.EPIC).

For estimateCellCounts2, I explain this here if you are using a 450K array GitHub FlowSorted.Blood.EPIC. If you are using EPIC you can read the help file of estimateCellCounts2. Although we believe our method improves the cell estimates, you can load the library and use minfi's "auto" selection or the same option on estimateCellCounts2 (that is your decision as a user of the library).

```
library(FlowSorted.Blood.EPIC)
FlowSorted.Blood.EPIC <-
    libraryDataGet("FlowSorted.Blood.EPIC")
FlowSorted.Blood.EPIC

```

Then you can run estimateCellCounts according to minfi's help.

Good luck.
> 2 votes

---
Tags: r
---


# Question
Title: WGCNA error: Length of colors vector not compatible with number of objects in 'order'
```
plotDendroAndColors(bwnet$dendrograms[[1]], cbind(bwnet$unmergedColors, bwnet$colors),
                    c("unmerged", "merged"),
                    dendroLabels = FALSE,
                    addGuide = TRUE,
                    hang= 0.03,
                    guideHang = 0.05)

```

I am getting the following error. colrs modules are not coming in dendrogram graph. I am using 4gb ram laptop. is this problem because of laptop or my mistake in understanding code.

```
Error in .plotOrderedColorSubplot(order = order, colors = colors,
rowLabels = rowLabels, : Length of colors vector not compatible with
number of objects in 'order'. madam getting this error.

```

My total genes are 15456. Under blockwisemodule I chose 5000

# Answer
> 1 votes
Please follow the code in WGCNA tutorial section 2c. Since the analysis ran in blockwise manner, you need to plot the dendrograms for each block separately and subset the module colors to each block. The example code is in the tutorial.

---
Tags: r, bioconductor, wgcna
---


title: What alternatives are there to Bracken?

text:
# Question
Title: What alternatives are there to Bracken?
Bracken is the recommended post-processer of Kraken2 outputs and specifically is expected to improve relative abundance estimates. However, Bracken's Github has over 100 outstanding issues and the code seems to rely on Python 2, which is outdated. While many of the github issues are likely not due to bugs, it does seem like the software does not seem to be supported strongly at this time. I am curious what other alternatives exist. Are folks moving from kraken2/bracken to another solution for quick taxonomic hypothesis generation of metagenomics data?

# Answer
No, Braken is the *tour de force* of first stage post-hoc Kraken2 analysis.

Whether the code is supported or not it corrects bias in Kraken2 population frequency distributions. I'm not aware of any theoretical concerns.

Generally, Python2 isn't a good idea, but it is overshadowed with the bigger concern of population bias, which is what Braken is there to adjust.
> 1 votes

---
Tags: metagenome
---


# Question
Title: A bash script for running on a bunch of bam files
I have some bam files in this directory

```
/data/Continuum/WES/results/

```

I want to run `GATK` mutation calling over bam files

I googled and I realised for one function, I can do this

```
cd /data/Continuum/WES/vcf/

for file in *.bam ; do ./gatk CollectSequencingArtifactMetrics -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta done;

```

I have five functions to run but I do not know how to include all in the script

These are my functions

```
./gatk CollectSequencingArtifactMetrics -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

./gatk GetPileupSummaries -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.table -V af-only-gnomad.hg38.vcf.gz -L resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

./gatk CalculateContamination -I NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.table -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.contamination.table

./gatk GetSampleName -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.sample_name

./gatk Mutect2 -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.vcf -tumor NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.sample_name --af-of-alleles-not-in-resource 2.5e-06 --germline-resource af-only-gnomad.hg38.vcf.gz -pon 1000g_pon.hg38.vcf.gz

```

I thought to use `|` to include five functions but I do not know how to change input and output for each iteration

I tried this for the first function,

```
    FILES="/data/Continuum/WES/testAlignmentBROADGenome2/results/.MarkedDup.bam*"
for f in $FILES
do
  ./gatk CollectSequencingArtifactMetrics -I $f -O $f --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta
done

```

but returned this error

```
SAMException: Cannot read non-existent file: file:///data/Continuum/WES/testAlignmentBROADGenome2/results/.MarkedDup.bam*

```

Anyone, please help me in develop this?

# Answer
There are a couple ways to approach this, one is to create output files based on the `BAM` input file name - the other approach is to make a new folder (based on the input `BAM` filename containing result files with the same name for each input. I generally use the first method, but either work.

Here is a script that may run the 5 `GATK` functions. It calculates a prefix to be used for derivative output files `BASE`, and should run your pipeline on all the `BAM` files within the `INPUT_FOLDER`.

```
#!/bin/bash

set -e

INPUT_FOLDER="/data/Continuum/WES/results/"

for BAM_FILE in "$INPUT_FOLDER/"*.bam
do
        BASE=$(basename "$BAM_FILE" .bam) # NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup

        # function 1
        ./gatk CollectSequencingArtifactMetrics \
                -I "$BAM_FILE" \
                -O "$BASE" \
                --FILE_EXTENSION .txt \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

        # function 2
        ./gatk GetPileupSummaries \
                -I "$BAM_FILE" \
                -O "$BASE.targeted_sequencing.table" \
                -V af-only-gnomad.hg38.vcf.gz \
                -L resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

        # function 3
        ./gatk CalculateContamination \
                -I "$BASE.targeted_sequencing.table" \
                -O "$BASE.targeted_sequencing.contamination.table" \

        # function 4
        ./gatk GetSampleName \
                -I "$BAM_FILE" \
                -O "$BASE.targeted_sequencing.sample_name"

        # function 5
        ./gatk Mutect2 \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta \
                -I "$BAM_FILE" \
                -O "$BASE".vcf \
                -tumor "$BASE.targeted_sequencing.sample_name" \
                --af-of-alleles-not-in-resource 2.5e-06 \
                --germline-resource af-only-gnomad.hg38.vcf.gz \
                -pon 1000g_pon.hg38.vcf.gz
done

```

It may need some adjustments - I did not test this, and used the `GATK` incantations you provided (thank you for including all of your steps!)

It includes the addition of the command `set -e`, which depending on GATKs return status may stop the script if any command fails - otherwise the script will happily run through a failed command.

You also may need to adjust the location of your output, and the directory the script is in upon execution, as you are running `GATK` from the current directory (`./gatk`)

In order to make these long lines readable, I use the method of placing a `\` at the end of a line - this allows the command to wrap to the next line. Be sure that the `\` character is at the end of your line for this to function properly.
> 6 votes

---
Tags: bash, gatk, wes, mutect2
---


# Question
Title: Feature-selection microbiome data to build a model to predict cancer
I want to build a machine learning model to predict colorectal cancer based on 16S rRNA microbiome data (stool samples).

I have filtered the data using filtering approach by Duvallet(removing samples with fewer than 100 reads and OTUs with fewer than 10 reads,as well as OTUs which were present in fewer than 1% within the study) What do I have to do after the preprocessing step? How can I do the feature selection for the phenotype data and the microbiome data?

# Answer
> 0 votes
I think the question needs improving. As the question stands pre-processing to feature selection ML (that's whats being asked).

* put the species diversity into a pandas dataframe the last column is the target (cancer vs healthy).
* assess a variety of ML algorithms (because they give different answers). Naive bayes, linear regression and K-nearest neighbour (KNN) are simple algorithms, easy to run and at this level fine.
* each have a "feature selection" - all features will be positive.

At a personal level I like KNN, but avoid LR and NB.

I personally would simply suggest unsupervised learning.

---
Tags: python, metagenome, machine-learning, microbiome
---


# Question
Title: Two variants associated with same chromosomal position in genotype data for one individual
I'm trying to play around with my family's raw genetic data from Ancestry.com so that I can do some genetic analyses, and I noticed in the AncestryDNA raw data files that some regions of the genome have more than one variant associated with them and I would like to know why this is. The end goal is get this into a regularly formatted VCF, but I feel that these double variants will probably mess things up. So for example one might see something similar to this in the tab-delimited file:

```
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

```

(*Header line*: rsid chromosome position allele1 allele2)

I understand that when rsid numbers were created, sometimes multiple rsid numbers were accidentally assigned to the same region of the genome, but I do not understand how two separate variants can be assigned to one single region in the genome. I could be wrong but as I understand it, I should only see one line here (either a "T T" or an "I I") because genotypes are reported in comparison to the human reference genome and so each position should uniquely refer to a single nucleotide in the said human reference genome.

But in this example case one line says there is a "T T" genotype at this position, and the other line says there is an insertion(see "How do I read my raw DNA data?" here). One possible explanation I thought of is if this would just be saying that a T was inserted, but if this is the case, then as I understand it the reference allele present in the individual should still reported, just in conjunction with the allele at the index. (So assuming an extra "A" was inserted after the T, which could be the reference or alternative allele, something like this should be recorded: rs6010164 22 50515270 TA TA).

Can anyone help me understand why there are two separate variants associated with chromosomal position? Thanks!

# Answer
> 2 votes
First of all, **take these data with a very large pinch of salt**. This sort of targeted analysis is not designed to produce high quality genetic data but, to give an idea of a sample's ancestry. Given a file like the one you describe with no quality information, you have no way of knowing how reliable any of the variants called actually is. From personal experience with another, similar company, I expect that they aren't reliable at all.

I have tested four individuals using a commercial tool like the one you describe and while I don't doubt the validity of their main results (those displayed on their website), their raw data showed combinations of variants that were extremely unlikely to be real given that the individuals were still breathing.

This isn't to say that such companies don't do their job well, only that their job is limited to what they actually display. When you look into their raw data, you're on your own. Presumably, they apply various filters before displaying the main results, but these are most likely not applied to the raw data. So if you want to do your own analyses, this is unlikely to be a reliable dataset. Also, I would strongly urge you to seek the help of an expert when interpreting any data you do have.

---

That said, the rsIDs are simply identifiers of a specific genetic variation, not of a location. Any genomic position can have many (essentially infinite) possible variations. For example, a C residue at position 100 of chromosome 1 could become a G, a T or an A. Each of these will have its own, separate rsID. Now, if you also consider insertions and deletions, you can see why you can have essentially limitless possible variations at that specific position. So no, having more than one rsID describing the same position is not strange at all.

The way these tests work is they look at a list of possible variants, each of which has its own rsID, and they report what they find for that rsID in the sample. Then, they will report what they find at the rsID's position for each allele. This is explained in the FAQ article you linked to:

> The possible observations are A for adenine, C for cytosine, G for guanine, T for thymine, I for insertion’ and ‘D for deletion, or 0 for missing data. Column one provides the identifier (including the #rsID where possible). Columns two and three contain the chromosome and basepair position of the variant using human reference 37.1 coordinates. Columns four and five contain the two alleles observed at this variant (genotype)

Let's have a look at the two rsIds you mention.

1. rs6010164: This is a single nucleotide variant (SNV) which describes a C residue where the reference genome has a T.
2. rs267607236: this is a deletion at position 50076841, so a very different variant to the above. In this one, the sequence TGAG is deleted.

Your results are:

```
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

```

So, your sample doesn't have rs6010164 since it has the same sequence at that position as the reference. However, it does have the insertion as indicated by the two Is.

The Is and Ds of the insertions and deletions (indels) are going to be the easy bit here. The only way you can parse this sort of file and convert it to VCF is to get the list of target SNV rsIDs, check what they describe as the reference and as the variant and then, for every one, check which of the two is shown in your file.

If you're not an expert in this field, you will probably be better off hiring someone to do this for you. There are actually companies that offer this sort of service and can take these tab-delimited files as input. I even work for one!

# Answer
> 2 votes
It's pretty common for SNPs and indels that map to the same position in the genome to be assigned different rsIDs, particularly when the indel spans multiple base pairs in the reference. You can get more info about these particular variants here and here. Or you can see how both of these variants relate to the reference sequence here.

Without info about the reference allele it's pretty confusing (particularly because they appear to be reporting the opposite strand from the database (G vs. C/T), but I would interpret this to mean that both alleles match the reference for rs267607236 (they have the insertion), but do not match the reference for rs6010164 (G-C pair rather than T-A). This is also the most common genotype at this position, being found in almost 80% of people in most populations.

# Answer
> 2 votes
You need to understand a few key points for your example given:

(a) Ancestry (and 23andMe) will use the letters I **and** D to indicate it is an InDel/DELINS. You have to read the rsID definition to determine if the variant is an Insert or Delete from the reference. And therefore which letter is used to mean positive / variant and which for negative / reference. They simply use the other letter to indicate negative / reference.

rs267607236 is defined as a deletion. So the I's mean the reference was detected in this case. Therefore the Ancestry file indicates no deletion of the 4 base-pairs at coordinates 50515271-4.

(b) Ancestry and dbSNP VCF annotation files uses a "left-aligned" coordinate for ALL InDel's and generally designate the coordinate BEFORE the variant.

This makes sense for Inserts. It gives the coordinate of the last base-pair in the reference just before the insert starts. That coordinate plus 1 is the next base-pair after the insert in the reference.

This same method is used for deletions. But because the deleted values have coordinates in the reference, this makes deletions appear "one off". You could use the starting coordinate of the deletion when defining the variant (e.g. `[TGAG/-]`). Otherwise, you need the coordinate before if defining the deletion as `[CTGAG/C]`.

rs267607236 is listed in the OP's information as coordinate 50515270. Because it is a deletion, the deleted sequence actually starts with the base-pair at coordinate 50515270+1 or 50515271.

(c) As it ends up, for rs267607236 in question here, the dbSNP VCF and the dbSNP online website entry are different (both latest release 156 in this case). This could happen but still be consistent because there is a `T` flanking the 4 base-pair deletion. The starting reference sequence `...C TGAGT A...` becomes `CTA` after the deletion. But which `T` remains is not really known. The VCF file defines the deletion by retaining the trailing `T`. The online website retains the leading `T`.

And in fact, it appears the website entry starting coordinate is incorrect. Given the way they defined this variant within the flanking sequence, the online website should use the coordinate 22:50515270 (like in the Ancestry file) and not 22:50515271 as it does now. This may have existed with that coordinate in an earlier release used by Illumina / Ancestry when defining the pipeline. But why the VCF and online website are different now is confusing.

There are numerous other nuances, conventions and actual errors that have to be understand to properly convert a microarray raw data file to a VCF.

Here are all the entries that we referenced / consulted brought together for your perusal:

```
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

GRCh37 Reference sequence around rs267607236:
TTAC TGAGT AAGA  (70-74)

dbSNP Release 156 GRCh37 VCF:
NC_000022.10    50515270        rs6010164       T       C
NC_000022.10    50515269        rs267607236     CTGAG   C       # Note: leading T deleted

dbSNP web page definition:
Alleles:T>C             for rs6010164)
TTAC [T/C] GAGTAAGAT
22:50515270 (GRCh37)
NC_000022.10:50515270:T:C

Alleles:TGAGT>T         for rs267607236
TTACT [GAGT/-] AAGAT                                            # Note: trailing T deleted
22:50515271 (GRCh37)                # for consistency with flanking sequence, this should be ..270
NC_000022.10:50515271:TGAGT:T

```

It is not clear to me what Ancestry does when the same rsID has an insert and a deletion at the same position. For example, the line from dbSNP VCF for rs768019142 is:
`NC_000001.11 10055 rs768019142 TA T,TAA`
which says there are two variants in the rsID. (This is common, especially for SNPs with multiple variants possible). A deletion at 10056 OR a single letter insert after 10056. To represent this, instead of a single letter for the reference, they have given two base-pairs as the reference sequence. They probably should have defined two rsIDs at two adjacent locations. UCSC's versions of the VCF files are filtered and cleanup to simplify items like this.

# Answer
> 1 votes
For different variation type at the same location, each variation type is being assigned its own rsID by the dbSNP maintainers.

In your example, one chromosome in the pair has a single-nucleotide variation (T allele), while the other chromosome has an insertion variation.

Quote:

> dbSNP only merges rs numbers that have an identical set of mappings to a genome and the same allele type (e.g. both must be the same variation type and share one allele in common). We therefore would not merge a SNP and an indel (insertion/deletion) into a single rs number (different variation classes) since they represent two different types of mutational "events".

In your case, rs6010164 is of the Single Nucleotide Variation (SNV) variation type, while rs267607236 is of the Small deletions or insertions (DELINS) variation type.

Quote:

> Despite its name, RefSNP is assigned to all variation types listed below with precise locations for both common and rare variations, including mutations. Most are typically small variations (\<= 50bp).
>
> * Single nucleotide variation (SNV)
> * Short multi-nucleotide changes (MNV)
> * Small deletions or insertions
> * Small STR repeats
> * retrotransposable element insertions

---
Tags: format-conversion, reference-genome, genotyping
---


# Question
Title: How to generate rsid values locally from VCF file
I recently used Dante to obtain a WGS and have downloaded several files produced during this process. There is a file that ends '.vcf' which I did think was a VCF file but I'm not sure now after a previous question. Part of that file is included below. Its about 1.1GB in size.

It appears to me to follow VCF format at least to a degree but it is missing the rsid's. The field where an rsid would usually be in a VCF file is replaced with a period (.).

I would like to backfill (annotate?) this file with known rsid's (where defined). There are nearly 4 million of them and online tools don't seem ideal, I would prefer to be able to do it locally, if that is possible. I am not looking for software to do this, but a data file which specifies which rsid is assigned to a particular chr:pos:allele. That may or may not be possible - perhaps I misunderstand and its far more complex than that. But from what I can tell, for a particular reference genome, and a particular build of dbsnp, there is a simple relationship as described.

I want to do this, for instance, to be able to use various tools, like SNPedia. Dante appears to use GRCh37, vs SNPedia using mainly GRCh38, which I can see is an issue, and will probably require me to buy a GRCh38 conversion from Dante (or do it myself). But please ignore that particular issue right now.

I don't quite understand what the dependencies are. I've seen comments that suggest that reference genome is not important to generating rsid's from VCF files and that only dbsnp build number needs to match the environment that the VCF file was generated in, but that does not seem to make much sense to me. From what I can see, if you are converting chr:pos information to rsid's, you need to be using the right reference genome. Does dbsnp build really matter? Should you not use the latest build, if it is continually being added to - surely rsid's are not reused? And then a post on the ANNOVAR site seems to suggest that no one really knows what rsid's are and that to different people they are different things, further muddying the waters.

So in terms of question - is there a file, downloadable from somewhere, that contains mappings from chr:pos to rsid, and if so, how does one select the correct file for a particular dataset (IE generated with a particular reference genome, with a particular dbsnp, to be used with references that have used a particular reference/dbsnp build like GRCh38 / dbsnp build 151)?

This file (header first), downloaded from Dante, is named .filtered.snp.vcf

```
    ##fileformat=VCFv4.2
    ##FILTER=<ID=PASS,Description="All filters passed">
    ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths (counting only informative reads out of the total reads) for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions for alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
    ##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
    ##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
    ##FORMAT=<ID=GP,Number=G,Type=Float,Description="Phred-scaled posterior probabilities for genotypes as defined in the VCF specification">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=MB,Number=4,Type=Integer,Description="Per-sample component statistics to detect mate bias">
    ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
    ##FORMAT=<ID=PRI,Number=G,Type=Float,Description="Phred-scaled prior probabilities for genotypes">
    ##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
    ##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias">
    ##FORMAT=<ID=SQ,Number=A,Type=Float,Description="Somatic quality">
    ##DRAGENCommandLine=<ID=HashTableBuild,Version="SW: 01.003.044.3.9.5, HashTableVersion: 8",CommandLineOptions="dragen --build-hash-table true --output-directory grch37 --ht-reference grch37.fasta">
    ##DRAGENCommandLine=<ID=dragen,Version="SW: 05.121.645.4.0.3, HW: 05.121.645",Date="Thu Mar 02 09:58:00 UTC 2023",CommandLineOptions="-f -c output.cfg -r /references/grch37">
    ##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
    ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
    ##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
    ##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
    ##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
    ##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
    ##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (informative and non-informative); some reads may have been filtered based on mapq etc.">
    ##INFO=<ID=FractionInformativeReads,Number=1,Type=Float,Description="The fraction of informative reads out of the total reads">
    ##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
    ##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
    ##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
    ##FILTER=<ID=DRAGENSnpHardQUAL,Description="Set if true:QUAL < 3">
    ##FILTER=<ID=DRAGENIndelHardQUAL,Description="Set if true:QUAL < 3">
    ##FILTER=<ID=LowDepth,Description="Set if true:DP <= 1">
    ##FILTER=<ID=PloidyConflict,Description="Genotype call from variant caller not consistent with chromosome ploidy">
    ##FILTER=<ID=base_quality,Description="Site filtered because median base quality of alt reads at this locus does not meet threshold">
    ##FILTER=<ID=filtered_reads,Description="Site filtered because too large a fraction of reads have been filtered out">
    ##FILTER=<ID=fragment_length,Description="Site filtered because absolute difference between the median fragment length of alt reads and median fragment length of ref reads at this locus exceeds threshold">
    ##FILTER=<ID=low_af,Description="Allele frequency does not meet threshold">
    ##FILTER=<ID=low_depth,Description="Site filtered because the read depth is too low">
    ##FILTER=<ID=low_frac_info_reads,Description="Site filtered because the fraction of informative reads is below threshold">
    ##FILTER=<ID=low_normal_depth,Description="Site filtered because the normal sample read depth is too low">
    ##FILTER=<ID=long_indel,Description="Site filtered because the indel length is too long">
    ##FILTER=<ID=mapping_quality,Description="Site filtered because median mapping quality of alt reads at this locus does not meet threshold">
    ##FILTER=<ID=multiallelic,Description="Site filtered because more than two alt alleles pass tumor LOD">
    ##FILTER=<ID=non_homref_normal,Description="Site filtered because the normal sample genotype is not homozygous reference">
    ##FILTER=<ID=no_reliable_supporting_read,Description="Site filtered because no reliable supporting somatic read exists">
    ##FILTER=<ID=panel_of_normals,Description="Seen in at least one sample in the panel of normals vcf">
    ##FILTER=<ID=read_position,Description="Site filtered because median of distances between start/end of read and this locus is below threshold">
    ##FILTER=<ID=RMxNRepeatRegion,Description="Site filtered because all or part of the variant allele is a repeat of the reference">
    ##FILTER=<ID=str_contraction,Description="Site filtered due to suspected PCR error where the alt allele is one repeat unit less than the reference">
    ##FILTER=<ID=too_few_supporting_reads,Description="Site filtered because there are too few supporting reads in the tumor sample">
    ##FILTER=<ID=weak_evidence,Description="Somatic variant score does not meet threshold">
    ##ML=PackageVersion=4,ModelVersion=4,ModelCRC=2921457924
    ##contig=<ID=1,length=249250621>
    ##contig=<ID=2,length=243199373>
    ##contig=<ID=3,length=198022430>
    ##contig=<ID=4,length=191154276>
    ##contig=<ID=5,length=180915260>
    ##contig=<ID=6,length=171115067>
    ##contig=<ID=7,length=159138663>
    ##contig=<ID=8,length=146364022>
    ##contig=<ID=9,length=141213431>
    ##contig=<ID=10,length=135534747>
    ##contig=<ID=11,length=135006516>
    ##contig=<ID=12,length=133851895>
    ##contig=<ID=13,length=115169878>
    ##contig=<ID=14,length=107349540>
    ##contig=<ID=15,length=102531392>
    ##contig=<ID=16,length=90354753>
    ##contig=<ID=17,length=81195210>
    ##contig=<ID=18,length=78077248>
    ##contig=<ID=19,length=59128983>
    ##contig=<ID=20,length=63025520>
    ##contig=<ID=21,length=48129895>
    ##contig=<ID=22,length=51304566>
    ##contig=<ID=X,length=155270560>
    ##contig=<ID=Y,length=59373566>
    ##contig=<ID=MT,length=16569>
    ##contig=<ID=GL000207.1,length=4262>
    ##contig=<ID=GL000226.1,length=15008>
    ##contig=<ID=GL000229.1,length=19913>
    ##contig=<ID=GL000231.1,length=27386>
    ##contig=<ID=GL000210.1,length=27682>
    ##contig=<ID=GL000239.1,length=33824>
    ##contig=<ID=GL000235.1,length=34474>
    ##contig=<ID=GL000201.1,length=36148>
    ##contig=<ID=GL000247.1,length=36422>
    ##contig=<ID=GL000245.1,length=36651>
    ##contig=<ID=GL000197.1,length=37175>
    ##contig=<ID=GL000203.1,length=37498>
    ##contig=<ID=GL000246.1,length=38154>
    ##contig=<ID=GL000249.1,length=38502>
    ##contig=<ID=GL000196.1,length=38914>
    ##contig=<ID=GL000248.1,length=39786>
    ##contig=<ID=GL000244.1,length=39929>
    ##contig=<ID=GL000238.1,length=39939>
    ##contig=<ID=GL000202.1,length=40103>
    ##contig=<ID=GL000234.1,length=40531>
    ##contig=<ID=GL000232.1,length=40652>
    ##contig=<ID=GL000206.1,length=41001>
    ##contig=<ID=GL000240.1,length=41933>
    ##contig=<ID=GL000236.1,length=41934>
    ##contig=<ID=GL000241.1,length=42152>
    ##contig=<ID=GL000243.1,length=43341>
    ##contig=<ID=GL000242.1,length=43523>
    ##contig=<ID=GL000230.1,length=43691>
    ##contig=<ID=GL000237.1,length=45867>
    ##contig=<ID=GL000233.1,length=45941>
    ##contig=<ID=GL000204.1,length=81310>
    ##contig=<ID=GL000198.1,length=90085>
    ##contig=<ID=GL000208.1,length=92689>
    ##contig=<ID=GL000191.1,length=106433>
    ##contig=<ID=GL000227.1,length=128374>
    ##contig=<ID=GL000228.1,length=129120>
    ##contig=<ID=GL000214.1,length=137718>
    ##contig=<ID=GL000221.1,length=155397>
    ##contig=<ID=GL000209.1,length=159169>
    ##contig=<ID=GL000218.1,length=161147>
    ##contig=<ID=GL000220.1,length=161802>
    ##contig=<ID=GL000213.1,length=164239>
    ##contig=<ID=GL000211.1,length=166566>
    ##contig=<ID=GL000199.1,length=169874>
    ##contig=<ID=GL000217.1,length=172149>
    ##contig=<ID=GL000216.1,length=172294>
    ##contig=<ID=GL000215.1,length=172545>
    ##contig=<ID=GL000205.1,length=174588>
    ##contig=<ID=GL000219.1,length=179198>
    ##contig=<ID=GL000224.1,length=179693>
    ##contig=<ID=GL000223.1,length=180455>
    ##contig=<ID=GL000195.1,length=182896>
    ##contig=<ID=GL000212.1,length=186858>
    ##contig=<ID=GL000222.1,length=186861>
    ##contig=<ID=GL000200.1,length=187035>
    ##contig=<ID=GL000193.1,length=189789>
    ##contig=<ID=GL000194.1,length=191469>
    ##contig=<ID=GL000225.1,length=211173>
    ##contig=<ID=GL000192.1,length=547496>
    ##contig=<ID=NC_007605,length=171823>
    ##contig=<ID=hs37d5,length=35477943>
    ##reference=file:///references/grch37/reference.bin
    ##bcftools_filterVersion=1.9+htslib-1.9
    ##bcftools_filterCommand=filter -i 'TYPE="snp" ' /processing/sample.raw.vcf.gz; Date=Thu Mar  2 10:28:45 2023
    #CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  GFX0xxxxxx
    1   12783   .   G   A   6.16    PASS    AC=1;AF=0.5;AN=2;DP=37;FS=0;MQ=11.6;MQRankSum=-2.137;QD=0.45;ReadPosRankSum=2.878;SOR=1.949;FractionInformativeReads=0.973  GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:10,26:0.722:36:4,14:6,12:4:38,0,1:6.1614,2.5231,7.0198:0,34.77,37.77:10,0,26,0:4,6,14,12
    1   13868   .   A   G   4.69    PASS    AC=1;AF=0.5;AN=2;DP=61;FS=2.275;MQ=16.39;MQRankSum=2.368;QD=0.36;ReadPosRankSum=-1.078;SOR=1.136;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:36,25:0.41:61:16,15:20,10:4:37,0,8:4.6918,2.1112,13.42:0,34.77,37.78:15,21,8,17:23,13,12,13
    1   13980   .   T   C   23.14   PASS    AC=2;AF=1;AN=2;DP=2;FS=0;MQ=5.7;MQRankSum=0;QD=6.5;ReadPosRankSum=0;SOR=2.303;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,2:1:2:0,1:0,1:19:61,24,0:23.143,20.887,0.056833:0,34.77,37.77:0,0,2,0:0,0,1,1
    1   14464   .   A   T   14.11   PASS    AC=1;AF=0.5;AN=2;DP=48;FS=7.243;MQ=22.38;MQRankSum=3.14;QD=0.86;ReadPosRankSum=2.527;SOR=0.267;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:11,37:0.771:48:8,16:3,21:7:77,6,0:14.113,0.98271,7.8591:0,34.77,37.77:3,8,19,18:4,7,20,17
    1   14907   .   A   G   3.67    PASS    AC=1;AF=0.5;AN=2;DP=77;FS=0.893;MQ=26.78;MQRankSum=-0.058;QD=0.41;ReadPosRankSum=5.956;SOR=0.532;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:48,29:0.377:77:31,13:17,16:4:36,0,22:3.6675,2.453,27.555:0,34.77,37.77:27,21,15,14:28,20,15,14
    1   14930   .   A   G   3.26    PASS    AC=1;AF=0.5;AN=2;DP=82;FS=3.051;MQ=27.04;MQRankSum=-0.363;QD=0.37;ReadPosRankSum=-0.277;SOR=0.754;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:51,31:0.378:82:30,16:21,15:3:35,0,21:3.2627,2.7897,26.606:0,34.77,37.77:28,23,14,17:29,22,17,14
    1   15211   .   T   G   6.66    PASS    AC=1;AF=0.5;AN=2;DP=70;FS=0;MQ=21.21;MQRankSum=0.533;QD=0.49;ReadPosRankSum=-0.905;SOR=0.743;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:28,42:0.6:70:12,22:16,20:6:40,0,11:6.6596,1.2235,15.27:0,34.77,37.77:15,13,23,19:15,13,17,25

```

Then a couple of chunks from random bits of the file:

offset 100000 lines:

```
    1   82279586    .   T   C   62.12   PASS    AC=2;AF=1;AN=2;DP=25;FS=0;MQ=250;MQRankSum=0;QD=4.89;ReadPosRankSum=0;SOR=0.941;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,25:1:25:0,10:0,15:52:100,55,0:62.124,52.173,2.8992e-05:0,34.77,37.77:0,0,11,14:0,0,10,15
    1   82280153    .   G   A   62.73   PASS    AC=2;AF=1;AN=2;DP=31;FS=0;MQ=250;MQRankSum=0;QD=4.53;ReadPosRankSum=0;SOR=1.044;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,31:1:31:0,15:0,16:57:101,61,0:62.731,57.663,9.319e-06:0,34.77,37.77:0,0,18,13:0,0,14,17
    1   82280683    .   C   T   54.83   PASS    AC=2;AF=1;AN=2;DP=17;FS=0;MQ=250;MQRankSum=0;QD=5.77;ReadPosRankSum=0;SOR=1.061;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,17:1:17:0,8:0,9:45:93,49,0:54.832,45.914,0.00012581:0,34.77,37.77:0,0,10,7:0,0,11,6
    1   82280870    .   T   C   41.24   PASS    AC=1;AF=0.5;AN=2;DP=19;FS=1.808;MQ=250;MQRankSum=3.592;QD=2.58;ReadPosRankSum=1.858;SOR=1.14;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:11,8:0.421:19:8,2:3,6:40:76,0,45:41.236,0.00039866,47.817:0,34.77,37.77:4,7,2,6:5,6,6,2
    1   82280975    .   A   T   40.52   PASS    AC=1;AF=0.5;AN=2;DP=24;FS=1.569;MQ=250;MQRankSum=4.07;QD=2.02;ReadPosRankSum=1.025;SOR=1.085;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,10:0.417:24:5,8:9,2:40:75,0,50:40.524,0.00040643,53.142:0,34.77,37.77:7,7,4,6:8,6,4,6
    1   82281492    .   T   C   44.14   PASS    AC=1;AF=0.5;AN=2;DP=22;FS=6.923;MQ=250;MQRankSum=3.873;QD=2.2;ReadPosRankSum=-0.334;SOR=2.064;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:13,9:0.409:22:9,4:4,5:44:79,0,51:44.143,0.00018328,54.3:0,34.77,37.77:6,7,7,2:6,7,6,3
    1   82282062    .   T   C   8.34    PASS    AC=1;AF=0.5;AN=2;DP=26;FS=1.974;MQ=247.84;MQRankSum=3.383;QD=1.04;ReadPosRankSum=-1.561;SOR=0.863;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:21,5:0.192:26:9,3:12,2:8:42,0,29:8.3401,0.69098,32.665:0,34.77,37.77:9,12,3,2:9,12,1,4
    1   82282378    .   C   T   42.61   PASS    AC=1;AF=0.5;AN=2;DP=24;FS=1.685;MQ=250;MQRankSum=4.07;QD=2.02;ReadPosRankSum=2.898;SOR=1.23;FractionInformativeReads=1  GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,10:0.417:24:8,6:6,4:43:77,0,57:42.606,0.0002423,60.191:0,34.77,37.77:8,6,7,3:10,4,7,3
    1   82282788    .   C   T   48.99   PASS    AC=1;AF=0.5;AN=2;DP=30;FS=0;MQ=245.71;MQRankSum=2.889;QD=1.67;ReadPosRankSum=1.975;SOR=0.914;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,16:0.533:30:6,9:8,7:46:84,0,46:48.985,0.00010665,49.217:0,34.77,37.77:8,6,10,6:6,8,7,9
    1   82282947    .   G   T   47.31   PASS    AC=1;AF=0.5;AN=2;DP=32;FS=10.606;MQ=250;MQRankSum=4.796;QD=1.56;ReadPosRankSum=0.68;SOR=0.725;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:15,17:0.531:32:6,8:9,9:46:82,0,48:47.312,0.00011286,51.287:0,34.77,37.77:11,4,7,10:8,7,8,9

```

offset 2000000 lines:

```
    8   142260995   .   T   C   65.76   PASS    AC=2;AF=1;AN=2;DP=28;FS=0;MQ=248.93;MQRankSum=0;QD=4.69;ReadPosRankSum=0;SOR=1.179;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,28:1:28:0,14:0,14:53:104,56,0:65.759,52.946,2.3297e-05:0,34.77,37.77:0,0,11,17:0,0,13,15
    8   142261389   .   T   C   66.03   PASS    AC=2;AF=1;AN=2;DP=28;FS=0;MQ=250;MQRankSum=0;QD=4.67;ReadPosRankSum=0;SOR=0.997;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,28:1:28:0,12:0,16:54:104,57,0:66.026,54.004,1.812e-05:0,34.77,37.77:0,0,16,12:0,0,14,14
    8   142261541   .   G   A   65.34   PASS    AC=2;AF=1;AN=2;DP=30;FS=0;MQ=250;MQRankSum=0;QD=4.51;ReadPosRankSum=0;SOR=1.765;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,30:1:30:0,14:0,16:52:103,55,0:65.339,52.429,2.5886e-05:0,34.77,37.77:0,0,9,21:0,0,17,13
    8   142261618   .   A   G   64.46   PASS    AC=2;AF=1;AN=2;DP=27;FS=0;MQ=250;MQRankSum=0;QD=4.57;ReadPosRankSum=0;SOR=1.302;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB:PS   1|1:0,27:1:27:0,14:0,13:53:102,56,0:64.465,53.434,2.1227e-05:0,34.77,37.78:0,0,10,17:0,0,15,12:142261618
    8   142261621   .   C   T   60.33   PASS    AC=2;AF=1;AN=2;DP=25;FS=0;MQ=250;MQRankSum=0;QD=4.7;ReadPosRankSum=0;SOR=1.051;FractionInformativeReads=0.96    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB:PS 1|1:0,24:1:24:0,13:0,11:48:98,51,0:60.326,48.044,7.2481e-05:0,34.77,37.78:0,0,10,14:0,0,13,11:142261618
    8   142261967   .   C   T   68.75   PASS    AC=2;AF=1;AN=2;DP=37;FS=0;MQ=250;MQRankSum=0;QD=4.26;ReadPosRankSum=0;SOR=0.859;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,37:1:37:0,17:0,20:56:107,60,0:68.745,56.755,9.8367e-06:0,34.77,37.77:0,0,20,17:0,0,18,19
    8   142262327   .   A   G   68.48   PASS    AC=2;AF=1;AN=2;DP=42;FS=0;MQ=250;MQRankSum=0;QD=4.13;ReadPosRankSum=0;SOR=0.788;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,42:1:42:0,21:0,21:56:106,60,0:68.476,56.786,9.8367e-06:0,34.77,37.77:0,0,22,20:0,0,11,31
    8   142262413   .   G   A   62.88   PASS    AC=2;AF=1;AN=2;DP=27;FS=0;MQ=231.86;MQRankSum=0;QD=4.75;ReadPosRankSum=0;SOR=1.302;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,27:1:27:0,16:0,11:51:101,54,0:62.881,51.183,3.5205e-05:0,34.77,37.77:0,0,10,17:0,0,9,18
    8   142263074   .   G   A   51.24   PASS    AC=2;AF=1;AN=2;DP=20;FS=0;MQ=143.27;MQRankSum=0;QD=5.31;ReadPosRankSum=0;SOR=0.892;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,20:1:20:0,9:0,11:36:89,39,0:51.237,35.702,0.0012015:0,34.77,37.77:0,0,11,9:0,0,10,10
    8   142263171   .   C   T   52.9    PASS    AC=2;AF=1;AN=2;DP=22;FS=0;MQ=122.2;MQRankSum=0;QD=4.81;ReadPosRankSum=0;SOR=0.693;FractionInformativeReads=0.909    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,20:1:20:0,11:0,9:35:91,38,0:52.899,34.894,0.0014297:0,34.77,37.77:0,0,10,10:0,0,9,11

```

# Answer
Download the dbSNP VCF file (latest version) and map to your input VCF using `bcftools annotate`. You'll want to match reference sequence version (hg19/hg38) and use the latest version of dbSNP to get the most up to date information.

dbSNP VCF FTP folder: https://ftp.ncbi.nih.gov/snp/latest\_release/VCF/
> 3 votes

# Answer
As a partial answer to the specific questions about definitions:

* rsIDs define a specific location and type of variation\*, but the location of that variation can change within different reference genomes. This may include a change in chromosome, as well as a change in position.
* Different dbSNP builds can add additional variation, but the same rsID will *always* refer to the same variant in *all* dbSNP builds that refer to that variation.
* rsIDs can be invalidated in new dbSNP builds for various reasons. It's a good idea to always use the latest dbSNP build when annotating variants.
* rsIDs are never recycled. There is a prescribed process for merging and splitting variants that ensures backwards compatibility. If a variation is merged, I think the approach is to always define the earliest-numbered variant as the primary variation. If a variation is split, I think the approach is to use a new rsID for each of the separate variations.

* by "variation", I'm referring to a type of modification in the DNA sequence, rather than a specific instance of that modification. This is explained in the dbSNP documentation (an "ss" is a submitted variant, a submission made to dbSNP):

> Each ss that maps to the same position and is the same variant type is assigned to an existing RefSNP (rs) or assigned a new one.

Although that isn't entirely clear, this suggests that if a new variant is submitted that falls at the same position of an exiting one and is of the same type as the existing one, then the new one may simply be added to the same rsID as the existing one. So, for example, if we have a single nucleotide G =\> C variant at position N, and a new G =\> T variant is reported at the same position, then the new one could have the same rsID as the old one.

---

In response to an additional comment / question:

> a post on the ANNOVAR site seems to suggest that no one really knows what rsid's are and that to different people they are different things, further muddying the waters.... Is this guy an authority on this or not?

I can't comment on authority for that guy; the authorship is unclear to me. What matters is how the dbSNP people understand dbSNP (i.e. the documentation link above), because they are the authority on it.

My understanding of it was helped because many years ago I attended a close-knit workshop that included a talk by Mike Feolo, and I had been working on a human dataset with many oddities, so I was able to ask lots of questions about corner cases. Here's the text from a section on rsIDs that I added into the appendix of my thesis; things have probably changed since then, but I expect the broad strokes are still correct:

> Mike Feolo gave a talk at the \[Working With The HapMap\] Cambridge course on the dbSNP database. He talked about how new Single Nucleotide Polymorphisms (SNPs) were submitted, compared against currently existing SNPs, and merged if they were already in the database. The database RefSeq(rs) numbers are updated about twice a year, and whenever there is a new build of the human genome sequence. As a result of these builds, some rs numbers are merged to account for increased knowledge of mutations, a process that I was not aware of. This has an impact on my research, because I have received data from an Illumina 317k SNPchip, which I compare to the HapMap dataset. The rs numbers from the Illumina chip are for a specific dbSNP build, which may not be fully consistent with the most recent build (from which the HapMap data is retrieved). Mike has recommended that both the build number and the rs number are used when referring to a specific mutation. Sharma Buch pointed out that it was possible to ask Illumina for information about the merged rs numbers for their chips, which would allow datasets from different builds to be used together.
> 2 votes

# Answer
Generally, consumer WGS companies do not deliver annotated VCF's. Annotations include rsIDs, gene names or similar area designations, and others. But this can be done using the NIH dbSNP VCF files and the bcftools annotate command. The annotation files at UCSC are modified for their GenomeBrowser and are often filtered and cleaned up. While better and more consistent, they are usually not the latest release available from NIH.

Just make sure to match the build model and sequence (chromosome) naming convention between your annotation file and the original VCF. Dante usually delivers using the hs37d5 reference model. So use the .25 dbSNP VCFs files. But dbSNP VCFs use accession names whereas Dante's hs37d5 reference model uses GRCh numeric ones. So you will have to rename sequence names in one file or the other to make it work. Also, hopefully, the vendor has left aligned the InDels as that is what the dbSNP file has them defined as and what bcftools expects as well. Dante traditionally uses the DRAGEN pipeline with their Illumina sequencers; which does this. But they have recently switched to MGI DNB sequencers and the megaBOLT pipeline for customers outside the USA. We are not yet sure how they are processing InDels. Finally, Dante has historically delivered a separate SNP and InDel file. You will often want to merge them or annotate them both separately then merge. bcftools can do this.

The dbSNP VCF file has well over 1 billion rsIDs now defined. And many fields of additional information that can be annotated in. Often, if you just want to add a simple subset of the data, you can process the dbSNP VCF file to a simpler tab-separated (TSV, tab) file of entries. It can cut down the annotation file from 25 GB to 5 GB or less. See the bcftools annotation page for hints on using that tool to do this as well.

Personally, I find it easier to use online tools or services that have been tweaked. Especially if they have nice browsers. Sequencing dot com, usegalaxy.org and iobio.io are a few that come to mind. Not all allow the export of the annotated VCF.

This and similar topics of utilizing your results are heavily discussed in the consumerWGS Facebook group; if you also wish to participate there.
> 1 votes

---
Tags: vcf, rsid
---


# Question
Title: How to loop over certain files within a directory that belongs to certain samples
*This question was also asked on Biostars*

I have two samples where the sequencing was done on 3-lanes paired end reads. The two samples are A44943\_1 and A44944\_1.

I would like to generate a for in loop to merge the R1 files of sample A44943\_1 from the 3-lanes, and then do the same for R2 files of this sample, then I want to merge R1 files of sample A44944\_1 from the 3-lanes and then do the same for R2 files of this sample. So I thought about for in loop. My fastq files looks like:

```
R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L001_R1_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L001_R2_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L002_R1_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L002_R2_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H7YKTDRXY_CGCAACTA-GAATCCGA_L003_R1_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H7YKTDRXY_CGCAACTA-GAATCCGA_L003_R2_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L001_R1_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L001_R2_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L002_R1_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L002_R2_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H7YKTDRXY_CACAGACT-TGGTACAG_L003_R1_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H7YKTDRXY_CACAGACT-TGGTACAG_L003_R2_trimmed.fastq

```

Could you advise me on a for in loop to do this task?

Thanks

# Answer
You can certainly use a for loop for this if desired, but an easier method is to use the `cat` command - something like this:

```
cat R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L00[123]_R1_trimmed.fastq > R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_R1_trimmed.fastq
cat R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L00[123]_R2_trimmed.fastq > R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_R2_trimmed.fastq
cat R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L00[123]_R1_trimmed.fastq > R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_R1_trimmed.fastq
cat R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L00[123]_R2_trimmed.fastq > R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_R2_trimmed.fastq

```
> 5 votes

---
Tags: linux
---


# Question
Title: Try to understand better how to set out an effective sample size within the MetagenomSeq package in microbiome
Cross-posted on biostars: https://www.biostars.org/p/9572744/

I am trying to understand better what number to be used and how to put it in the function MRcoefs within the Metagenomeseq package in R for calculating differentially expressed OTUs between two condition: my code looks like:

```
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = 5.5)

```

Here I am setting the eff argument to 5.5 because I have calculated the effectiv\_size using this code:

```
Effective_size(fit)

```

which result in table of effective size for each feature (OTU), and then I got the mean of these

```
mean(Effective_size(fit))

```

which was 5.5 When I run this, it produced an error: prop outside \[0,1\], I understand it should be between 0-1, but why ? is not that number in eff is the mean that I calculated which bring 5.5 ? Also when I tried to run it with .5 as the document suggested as below :

```
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = .5)

```

it run nicely, no error , but it produced the same OTUs as when I did not include eff, I was expecting that some OTUs will be removed as they were below this threshold

Could you please help me ?

# Answer
> 0 votes
Set `numberEff = True`

```
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = 5.5, numberEff = True)

```

I suspect you've got the quantile value set, hence it will throw an error above 1.

---
Tags: r, microbial-genomics, 16rrna, microbiome
---


# Question
Title: Paired end reads with different read lengths
I downloaded some raw RNA-seq reads from NCBI using the sra-toolkit (Accession numbers: SRR2053159-64). The reads were generated using ABI SOLiD, and in colour space. I downloaded them in base space in fastq format using the `-B` option in `fastq-dump`. But after that, when I am looking at the data after running FastQC, I find that the forward read length is 50 while the reverse reads length is 35. How is this possible? This is causing all the reverse reads to fail the filtering criteria set in Trimmomatic `MINLEN:36`. What is the correct way to work around this?

I have already set `LEADING:20 TRAILING:20 AVGQUAL:25 SLIDINGWINDOW:10:30 MINLEN:36` Would it be a good idea to entirely remove the `MINLEN:36`, or would it be better if I make the threshold 25?

Actually I am trying to compare the transcriptomes of different bacteria under a certain stress using raw reads from NCBI. While for the rest of the set of bacteria the sequencing was done using either Illumina or Solexa, this one was done using ABI SOLiD. Is it a good idea to include this ABI SOLiD data along with the others for analysis?

# Answer
> 2 votes
> I find that the forward read length is 50 while the reverse reads length is 35. How is this possible?

Forward cycle length was longer on SOLiD than reverse cycle length; this is common for SOLiD data. If you want to work with your data, you'll need to adjust your trimming settings appropriately to account for this (e.g. getting rid of the minimum length entirely, or setting it to a lower value). Trimmomatic does have other trimming options that are not length-based. Given the options that you've mentioned, it should be fine to just drop that minimum length parameter.

> Is it a good idea to include this ABI SOLiD data along with the others for analysis?

No, especially not if you're converting to base space, then mapping. SOLiD data should be mapped in colorspace, so that small changes in base transitions don't have a huge effect on the mapped sequence. SOLiD data converted to base space sequence is not appropriate for mapping because it mixes different error modes (i.e. sequencing error, and base transition error). Bowtie (as in Bowtie v1) can map to a colorspace reference, and will correct errors based on that reference in the resultant SAM output. This means that if you trust your reference, you can use Bowtie as a data cleaner for colorspace data.

However, on a practical level, the output of ABI SOLiD machines for existing datasets was low enough that (coupled with the bioinformatics issues) it makes more financial sense to resequence samples than to try to navigate around the pitfalls of colorspace representation.

I have previously offered the following steps for dealing with colorspace data:

1. Transfer all the colorspace files onto an external hard disk
2. Delete all other copies of the colorspace files
3. Remove the hard drive from the computer
4. Use a sledgehammer or similar to squash the disk platters closer together
5. Withdraw \\$500 from the bank
6. Place the \\$500 on top of the hard drive
7. Return the hard drive (with the money) back to the client
8. Report to the client that there was insufficient data for a suitable analysis, and recommend that the experiment is repeated using a base-space sequencer

## A More detailed Explanation

*\[from here\]*

In colorspace sequences that are converted to base space, any error in the sequence will cause all the following bases to be incorrect. Any colour-space to base-space conversion needs to take into account (and correct) errors so that the base-space sequences are correct. When there is a sequence difference, the conversion needs to make sure that only that position is changed in the base-space version.

Consider the following sequences that map to the same position:

```
reference: G101320112
sequence1: X101120112
sequence2: X101312011
sequence3: X201320312

```

\[I chucked an INDEL in there to make things a bit harder\]

A naive base-space conversion would convert these sequences as follows:

```
reference: GTTGCTTGTC
sequence1: GTTGTCCACT
sequence2: GTTGCAGGTG
sequence3: GAACGAATGA

```

\[apologies if my conversion is incorrect. Fixes appreciated\]

Very similar colour-space sequences, but very different base-space sequences.

A more correct conversion would notice where the errors were in the sequences relative to the index, and modify the next colour-space base as well to something that looks appropriate:

```
reference: G101320112
sequence1: X101100112
sequence2: X101313011
sequence3: X231320332

```

This would end up with these converted sequences:

```
reference: GTTGCTTGTC
sequence1: GTTGTTTGTC
sequence2: GTTGCATTGG
sequence3: GATGCTTATC

```

Which look considerably better.

I hate colorspace because the conversions are very unintuitive, and difficult to explain to other people. About the only nice thing is that reverse complement is just the reverse, but this also means that aligners and assemblers need to be modified to account for that when working in colorspace (or double-encoded colorspace), and you can get weird unexpected chimeras (e.g. poly-A tails and poly-T heads merging). You can save a lot of pain and confusion by sticking with a base-space sequencer.

---
Tags: rna-seq, trimming
---


# Question
Title: Pathway Enrichment Analysis for one specific pathway?
I have a gene set containing a set of 1500 genes in which I did differential gene expression using limma voom pipeline. After this step I wanted to analyze a specific set of pathways the metabolism set from KEGG database.

```
hs_kegg_df <- msigdbr(species = "Homo sapiens") %>%
 dplyr::filter(
gs_cat == "C2", # This is to filter only to the C2 curated gene sets
gs_subcat == "CP:KEGG" # This is because we only want KEGG pathways
)

```

Then since I only want to analyze metabolism I did this (keep only the pathways envolved in metabolism):

```
hs_kegg_df_with_gs_source <- subset(hs_kegg_df, grepl("^hsa00|^hsa01|^hsa02", gs_exact_source))

```

My problem is that, for example, for the pathway KEGG arginine and proline metabolism the set contains 54 genes in which my kit only covers 4, and these for 4 genes are Differentially expressed yet when doing: ( I decided to do ORA because with 1500 genes GSEA is not recommended).

```
em <- enricher(genes, minGSSize=1, maxGSSize = 20,  universe = names(gene_list), TERM2GENE=hs_kegg_df,pAdjustMethod = "fdr")

```

This set does not appear with significant p value. I also tried passing to the enricher just that pathway in particular, the p value is 1. Because in hypergeometric test, the background genes are the same as the as DGE genes in all pathways giving this way 1-0=1.

How can I analyze just one pathway in particular given this specific conditions?

These are the 4 genes DGE and included in my kit:

```
  gs_name                              gene_symbol
  <chr>                                <chr>
1 KEGG_ARGININE_AND_PROLINE_METABOLISM ALDH2
2 KEGG_ARGININE_AND_PROLINE_METABOLISM ARG2
3 KEGG_ARGININE_AND_PROLINE_METABOLISM NOS1
4 KEGG_ARGININE_AND_PROLINE_METABOLISM ODC1

```

And this are all the 59.

here all the genes dge and all the genes in pathway:

```
 genes_dge <- c("ALDH2", "ARG2", "NOS1", "ODC1")
 all_genes <- c("ACY1", "AGMAT", "ALDH18A1", "ALDH1B1", "ALDH2", "ALDH3A2", "ALDH4A1", "ALDH7A1", "ALDH9A1", "AMD1", "AOC1", "ARG1", "ARG2", "ASL", "ASS1", "AZIN2", "CKB", "CKM", "CKMT1A", "CKMT1B",s"CKMT2", "CPS1", "DAO", "GAMT", "GATM", "GLS", "GLS2", "GLUD1", "GLUD2", "GLUD2", "GLUL", "GOT1", "GOT2", "LAP3", "MAOA", "MAOB", "NAGS", "NOS1", "NOS2", "NOS3", "OAT", "ODC1", "OTC", "P4HA1", "P4HA2", "P4HA3", "PRODH", "PRODH2", "PYCR1", "PYCR2", "PYCR3", "PYCR3", "SAT1","SAT2", "SMS", "SRM")

```

# Answer
Two questions

1. Significance
2. Lack of consistency between enzymes (log2) within the arginine and proline metabolism pathway.

. 1. We've dealt with significance in a past post on Limma (it was not voom though). I'll try and find the post and link it here.

2. The answer is straight forward: whilst all four enzymes include Arginine and proline metabolism pathway

* only ARG and ODC1 share exactly the same pathways. Thus these should give very similar results.
* NOS1 is a key enzyme in 6 different pathways.
* ALDH2 is involved in 13 different pathways.

Therefore NOS1 and ALDH2 would have "confounders" (from an experimental perspective, not host metabolism) because they are more complex than singly arginine and proline metabolism pathways.

The other factor is in cancer patients the transcriptomics will get messed up, ALDH2 could go all over the place.

*Prediction* ARG = ORD \< NOS1 \< ALDH2

```
**ALDH2** aldehyde dehydrogenase 2
Pathways
Glycolysis / Gluconeogenesis
Ascorbate and aldarate metabolism
Fatty acid degradation
Valine, leucine and isoleucine degradation
Lysine degradation
**Arginine and proline metabolism**
Histidine metabolism
Tryptophan metabolism
beta-Alanine metabolism
Glycerolipid metabolism
Pyruvate metabolism
Pantothenate and CoA biosynthesis

```

```
**ARG** arginase
Pathways
**Arginine and proline metabolism**
Glutathione metabolism

```

```
**NOS1** nitric oxide synthase
Pathways
Arginine biosynthesis
**Arginine and proline metabolism**
Calcium signaling pathway
Phagosome
Apelin signaling pathway

```

```
**ODC1** ornithine decarboxylase 1
Pathways
**Arginine and proline metabolism**
Glutathione metabolism

```
> 1 votes

---
Tags: differential-expression, limma, gsea, clusterprofiler
---


title: low mapping percentage in bwa mem

text:
# Question
Title: low mapping percentage in bwa mem
After aligning paired-end reads to a reference genome, I am getting low percentage:

```
34636166 + 0 in total (QC-passed reads + QC-failed reads)
34346076 + 0 primary
0 + 0 secondary
290090 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
5706057 + 0 mapped (16.47% : N/A)
5415967 + 0 primary mapped (15.77% : N/A)
34346076 + 0 paired in sequencing
17173038 + 0 read1
17173038 + 0 read2
3807784 + 0 properly paired (11.09% : N/A)
3996748 + 0 with itself and mate mapped
1419219 + 0 singletons (4.13% : N/A)
177284 + 0 with mate mapped to a different chr
65388 + 0 with mate mapped to a different chr (mapQ>=5)

```

What has to be done to increase the mapping percentage?

I'm attaching the fastqc report before trimming and after trimming too first- R1 without trimming second- R2 without trimming

third- forward paired fourth- reverse paried

rest are overrepresentations

After using trimmomatic the quality is getting more bad and even the mapping is of low percentage

Mash screen result in the last image

# Answer
Just answering, following extensive back and forth in the comments.

I think that your reads are fine:

1. FASTQC doesn't see anything too weird. You see some noise at the first few basecalls but I am not too alarmed about that.
2. FASTQC finds overrepresented only a telomeric repeat and its reverse complement (TTAGG). That is something you'd expect to find in WGS reads.
3. Mash screen finds that the only match, with a small but non-trivial number of k-mer matches, is the monarch butterfly *Danaus plexippus*, e.g. another lepidoptera. So that makes sense.

**Therefore, the issue is likely that your reference genome has issues or has been swapped out with something else.**

Where did your reference assembly come from?

* Did it get mixed up somehow?
* What does it look like when you run QUAST?
* What does it look like when you run `mash dist` of your assembly against the same mash database? Is it maybe a bunch of endosymbionts or other contaminants?
* What happens if you map your reads against a high-quality lepidopteran assembly like the monarch butterfly or this?
* What happens if you map your assembly against a high-quality lepidopteran assembly? Is it syntenic?
> 1 votes

---
Tags: sam, samtools, reference-genome, bwa
---


# Question
Title: Why is there gene without annotated location?
I came across this gene **Pepc** (peptidase C)

NCBI gene database has a page for it here https://www.ncbi.nlm.nih.gov/gene/109616

But it does not contain any information on its coding sequence or genomic location. It just says it sits on chromosome 1.

Can someone familiar with how the reference genome is annotated explain why this is?

# Answer
This gene/protein appears to have been described exclusively as an isozyme / electrophoretic mobility assay of a protein.

See for example this PepC mutant page:

> The Pepc locus controls electrophoretic mobility of the enzyme peptidase C of red blood cells. This LT/J congenic strain carries the Pepca allele from C57BL/6 that determines a slow anodally migrating electrophoretic band.

The relevant publications describing mutants are from the 70s and 80s, e.g. here, when there were no genomic resources for mice.

However, linkage mapping of electrophoretic variants could be used to establish linkage to loci on chromosome 1, apparently.

Very likely PepC is encoded by *some* annotated gene on chr1, but no one knows which one exactly because they don't really know much about the protein other than its electrophoretic variability. But that's just a guess by me.

When the databases were set up, they had to accommodate 1) mouse mutant strains, and 2) published articles referencing the protein, so they did their best.

However, no one ever followed up to link everything together, so we just have these orphan entries.
> 1 votes

---
Tags: sequence-annotation
---


title: Analysis of RNA sequencing with zUMIs

text:
# Question
Title: Analysis of RNA sequencing with zUMIs
I wrote the .yaml file which was according the instruction from zUMIs:

```
###########################################
#Welcome to zUMIs
#below, please fill the mandatory inputs
#We expect full paths for all files.
###########################################

#define a project name that will be used to name output files
project: SmartSeq3

#Sequencing File Inputs:
#For each input file, make one list object & define path and barcode ranges
#base definition vocabulary: BC(n) UMI(n) cDNA(n).
#Barcode range definition needs to account for all ranges.You can  give several comma-separated ranges for BC & UMI sequences, eg.BC(1-6,20-26)
#you can specify between 1 and 4 input files

sequence_files:
  file1:
    name: /home/isidora/IRCCS_Candiolo/Fastq/*_R1_001_val_1.fq ##/home/isidora/IRCCS_Candiolo/Fastq/zUMIsFiles1/ #path to first file
    base_definition:
      - cDNA(1-150)
      - BC(1-8)
      - UMI(12-19)
  file2:
    name: /home/isidora/IRCCS_Candiolo/Fastq/*_R2_001_val_2.fq ## /home/isidora/IRCCS_Candiolo/Fastq/zUMIsFiles2/ #path to second file
    base_definition:
      - cDNA(1-150)
      - BC(1-8)
      - UMI(12-19)

#reference genome setup
reference:
  STAR_index: /home/isidora/IRCCS_Candiolo/Fastq/ #path to STAR genome index
  GTF_file: /home/isidora/IRCCS_Candiolo/Fastq/hg38.refGene.gtf #path to gene annotation file in GTF format
  exon_extension: no #extend exons by a certain width?
  extension_length: 0 #number of bp to extend exons by
  scaffold_length_min: 0 #minimal scaffold/chromosome length to consider (0 = all)
  additional_files: /home/isidora/IRCCS_Candiolo/Fastq/hg38.fa #Optional parameter. It is possible to give additonal reference sequences here, eg ERCC.fa
  additional_STAR_params: pGe.sjdbOverhang > 0 #Optional parameter. you may add custom mapping parameters to STAR here
  #output directory
out_dir: /home/isidora/IRCCS_Candiolo/Fastq/zUMIs/ #specify the full path to the output directory

###########################################
#below, you may optionally change default parameters
###########################################

#number of processors to use
num_threads: 20
mem_limit: null #Memory limit in Gigabytes, null meaning unlimited RAM usage
#barcode & UMI filtering options
#number of bases under the base quality cutoff that should be filtered out.
#Phred score base-cutoff for quality control.
filter_cutoffs:
  BC_filter:
    num_bases: 3
    phred: 20
  UMI_filter:
    num_bases: 2
    phred: 20

#Options for Barcode handling
#You can give either number of top barcodes to use or give an annotation of cell barcodes
#If you leave both barcode_num and barcode_file empty, zUMIs will perform automatic cell barcode selection for you!
barcodes:
  barcode_num: null
  barcode_file: /home/isidora/IRCCS_Candiolo/Fastq/SampleSheet.txt
  barcode_sharing: null #Optional for combining several barcode sequences per cell (see github wiki)
  automatic: yes #Give yes/no to this option. If the cell barcodes should be detected automatically. If the barcode file is given in combination with automatic barcode detection, the list of giiven barcodes will be used as whitelist.
  BarcodeBinning: 1 #Hamming distance binning of close cell barcode sequences.
  nReadsperCell: 100 #Keep only the cell barcodes with atleast n number of reads.
  demultiplex: yes #produce per-cell demultiplexed bam files.

#Options related to counting of reads towards expression profiles
counting_opts:
  introns: yes #can be set to no for exon-only counting.
  intronProb: no #perform an estimation of how likely intronic reads are to be derived from mRNA by comparing to intergenic counts.
  downsampling: '0' #Number of reads to downsample to. This value can be a fixed number of reads (e.g. 10000) or a desired range (e.g  10000-20000) Barcodes with less than <d> will not be reported. 0 means adaptive downsampling. Default: 0
  strand: 0 #Is the library stranded? 0 = unstranded, 1 = positively stranded, 2 = negatively stranded
  Ham_Dist: 1 #Hamming distance collapsing of UMI sequences.
  velocyto: no #Would you like velocyto to do counting of intron-exon spanning reads
  primaryHit: yes #Do you want to count the primary Hits of multimapping reads towards gene expression levels?
  multi_overlap: no #Do you want to assign reads overlapping to multiple features?
  fraction_overlap: 0 #minimum required fraction of the read overlapping with the gene for read assignment to genes
  twoPass: yes #perform basic STAR twoPass mapping

#produce stats files and plots?
make_stats: yes

#Start zUMIs from stage. Possible TEXT(Filtering, Mapping, Counting,Summarising). Default: Filtering.
which_Stage: Filtering

#define dependencies program paths
samtools_exec: samtools #samtools executable
Rscript_exec: Rscript #Rscript executable
STAR_exec: STAR #STAR executable
pigz_exec: pigz #pigz executable

#below, fqfilter will add a read_layout flag defining SE or PE

```

When I start the program with this code `zUMIs/zUMIs.sh -c -y ~/IRCCS_Candiolo/Fastq/zUMIs/zUMIs.yaml`, I get errors:

```
Using miniconda environment for zUMIs!
 note: internal executables will be used instead of those specified in the YAML file!

 You provided these parameters:
 YAML file:     /home/isidora/IRCCS_Candiolo/Fastq/zUMIs/zUMIs.yaml
 zUMIs directory:               /home/isidora/IRCCS_Candiolo/Fastq/zUMIs
 STAR executable                STAR
 samtools executable            samtools
 pigz executable                pigz
 Rscript executable             Rscript
 RAM limit:   null
 zUMIs version 2.9.7e

Wed Aug 23 12:06:25 CEST 2023
WARNING: The STAR version used for mapping is 2.7.3a and the STAR index was created using the version 2.7.4a. This may lead to an error while mapping. If you encounter any errors at the mapping stage, please make sure to create the STAR index using STAR 2.7.3a.
Filtering...
sh: 1: Syntax error: Unterminated quoted string
sh: 1: Syntax error: Unterminated quoted string
Wed Aug 23 12:06:26 CEST 2023
Error in eval(bysub, parent.frame(), parent.frame()) :
  object 'XC' not found
Calls: cellBC -> [ -> [.data.table -> eval -> eval
In addition: Warning message:
In data.table::fread(bccount_file, header = FALSE, col.names = c("XC",  :
  File '/home/isidora/IRCCS_Candiolo/Fastq/zUMIs//SmartSeq3.BCstats.txt' has size 0. Returning a NULL data.table.
Execution halted
Mapping...
[1] "2023-08-23 12:06:28 CEST"
Warning message:
In data.table::fread(cmd = paste(samtools, "view", filtered_bams[1],  :
  File '/tmp/RtmpKSSqU2/file57f5393a8fde' has size 0. Returning a NULL data.table.

EXITING because of fatal PARAMETERS error: pGe.sjdbOverhang <=0 while junctions are inserted on the fly with --sjdbFileChrStartEnd or/and --sjdbGTFfile
SOLUTION: specify pGe.sjdbOverhang>0, ideally readmateLength-1
Aug 23 12:07:18 ...... FATAL ERROR, exiting
Wed Aug 23 12:07:18 CEST 2023
Counting...
[1] "2023-08-23 12:07:30 CEST"
Error in fread(paste0(opt$out_dir, "/zUMIs_output/", opt$project, "kept_barcodes_binned.txt")) :
  File '/home/isidora/IRCCS_Candiolo/Fastq/zUMIs//zUMIs_output/SmartSeq3kept_barcodes_binned.txt' does not exist or is non-readable. getwd()=='/home/isidora/IRCCS_Candiolo/Fastq/zUMIs'
Execution halted
Wed Aug 23 12:07:30 CEST 2023
Loading required package: yaml
Loading required package: Matrix
[1] "loomR found"
Error in gzfile(file, "rb") : cannot open the connection
Calls: rds_to_loom -> readRDS -> gzfile
In addition: Warning message:
In gzfile(file, "rb") :
  cannot open compressed file '/home/isidora/IRCCS_Candiolo/Fastq/zUMIs//zUMIs_output/expression/SmartSeq3.dgecounts.rds', probable reason 'No such file or directory'
Execution halted
Wed Aug 23 12:07:32 CEST 2023
Descriptive statistics...
[1] "I am loading useful packages for plotting..."
[1] "2023-08-23 12:07:33 CEST"
Error in data.table::fread(paste0(opt$out_dir, "/zUMIs_output/", opt$project,  :
  File '/home/isidora/IRCCS_Candiolo/Fastq/zUMIs//zUMIs_output/SmartSeq3kept_barcodes.txt' does not exist or is non-readable. getwd()=='/home/isidora/IRCCS_Candiolo/Fastq/zUMIs'
Execution halted

```

I am new at bioinformatic field and I do not understand the meaning of these errors and the things I should correct. Can anyone help me please?

# Answer
I'd advise you to keep away from ZUMIs unless you're doing the same workflow that it was designed for. I spent many days trying to sort out issues with that program, and eventually gave up; it wasn't created with the expectation of errors or unexpected input, and had too many programming minefields for the data I was working with.

ZUMIs should be stopping on the first error and explaining what's wrong. As demonstrated by your situation, it's not; it's trying to keep on processing, barreling through and dragging errors along for the ride. It's not exactly beginner-friendly code.

In the end, I changed to using Salmon/Alevin, which had a more fail-resistant interface, and (eventually) similar barcode flexibility:

https://combine-lab.github.io/alevin-fry-tutorials/2021/custom-geo/
> 3 votes

---
Tags: command-line, bash
---


# Question
Title: GTF To Use With Broad Institute Fasta In STAR
To make a STAR reference genome used for alignment you need a genome fasta (the resources fasta provided by the Broad Institute in my case) and a GTF fasta. I tried to use the Ensembl human genome fasta fro hg38 but that does not work. I have looked in the GATK resources package but there is no `.gtf` file there. I really don't know what to do in this case? Is there anyone who can provide some guidance as to what `.gtf` file I should use?

# Answer
If you are looking for a genome and its corresponding annotations, you can use the following, it is from `nf-core` (link to the original doc):

```
#!/bin/bash

VERSION=108
wget -L ftp://ftp.ensembl.org/pub/release-$VERSION/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz
wget -L ftp://ftp.ensembl.org/pub/release-$VERSION/gtf/homo_sapiens/Homo_sapiens.GRCh38.$VERSION.gtf.gz

```

BTW the most up-to-date version as of today is 109.
> 1 votes

# Answer
I've had success in the past generating STAR indexes using the Gencode GTF files available from here:

https://www.gencodegenes.org/human/

* Comprehensive gene annotation (ALL regions)
* Genome sequence (ALL)

If you make sure that you're getting both the genome and the GTF files from the same source, it's more likely that they will work well together.

Note that the ALL annotation includes assembly patches and alternate haplotypes. The alternate haplotypes mean that population variation is captured, which usually increases the alignment rate of reads.

Unfortunately, most bioinformatics workflows don't properly handle alternate assemblies. As one example, reads aligning to those areas may end up classified as mapping to multiple genomic locations, even though they actually aren't. So it's a tradeoff of a more accurate alignment vs more hassle dealing with those alignments downstream.
> 1 votes

---
Tags: star
---


# Question
Title: Scale on the y-axis in a tree rotated with 'layout_dendrogram()' using ggtree
I'm editing a time-scaled tree using ggtree. The dates are correctly adjusted when I apply the 'mrsd' command and 'theme\_tree2()', but when I try to rotate the tree with 'layout\_dendrogram()', the scale still remains on the x-axis with no dates. Is there a way to place the date scale on the y-axis in a tree rotated with 'layout\_dendrogram()'?

```
# Code without layout_dendrogram()
ggtree(my_tree, mrsd='2022-12-31') %<+%
  metadato_arbol +
  geom_tippoint(size = 2) +
  theme_tree2()

```

```
# Using layout_dendrogram()
ggtree(my_tree, mrsd='2022-12-31') %<+%
  metadato_arbol +
  geom_tippoint(size = 2) +
  theme_tree2() + layout_dendrogram()

```

# Answer
Maybe try the following via gridExtra,

```
library(gridExtra)
library(ggtree)
library(ape)
set.seed(2023-08-23)
my_tree <- rtree(30)
p <- ggtree(my_tree, mrsd='2022-12-31') %<+%
  metadato_arbol +
  geom_tippoint(size = 2) +
  theme_tree2() # + geom_text(aes(label=node))
gridExtra::grid.arrange(p, p %>% rotate(33) %>% rotate(52), ncol=2)

```
> 0 votes

---
Tags: r, phylogenetics, ggtree
---


title: Model optimization causes overfitting

text:
# Question
Title: Model optimization causes overfitting
I'm working on an ML model to predict colorectal cancer using data from the human microbiome (collected from fecal samples).

My dataset is made up of metadata (containing phenotypic information on the N=156 population on which the study is being carried out,  and an OTU table containing a taxonomy of bacteria found and their abundance for each sample from one person (each person has one and only one sample).

I filtered the OTU table by eliminating sparse data (I eliminated OTUs below a certain threshold) after which I calculated the relative abundance for each OTU.
I deleted samples with missing values.
I made the metadata join (Subject ID Sample ID Age (years) Gender BMI (kg/m²) Country of Residence Diagnosis TNM Stage AJCC Stage Localization) and the OTU table.

In the end, I ended up with 114 rows × 221 columns .

I normalized the data and then built an SVM(SVC) model with the rbf kernel with C=0.05 and gamma=0.001 I got

```
Training Accuracy: 0.967032967032967
Test Accuracy: 0.782608695652174
SVC(AUC=0.64)

```

When I try to optimize the model using the code below, the values I get for gamma and C give me overfitting (trainset accuracy is 100%).

```
# import packages
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV
# range of values to try for gamma and C
Cs = [0.00001,0.0001,0.001,0.1,1,2,5,10,100,1000]
gammas = [0.00001,0.0001,0.001,0.1,1,10,100,10000]
# number of segments for cross validation
nfolds = 5 # number of segments for class validation
param_grid = {'C':Cs ,'gamma': gammas}
# search for the best parameters using cross-validation
grid_search = GridSearchCV(SVC(kernel='rbf'),param_grid,cv=nfolds)
grid_search.fit(X_train,y_train)
grid_search.best_params_

```

# Answer
> 0 votes
You are missing a validation split triggering data leakage. The error is here:

```
grid_search.fit(X_train,y_train)

```

You can't use `X_train` and `y_train` here. The correction is

```
grid_search.fit(X_val,y_val)

```

Generating `X_val` and `y_val` is described below.

The code used to generate

> Training Accuracy: 0.967032967032967 Test Accuracy: 0.782608695652174

... included one `test_split_train` and you need to perform a second `test_split_train` to obtain `X_val`, and `y_val`. The original code gave a pretty good as a result, the training score will always be higher.

The diagnosis of the code presented is there is 100% data leakage. What is missing is a second `test_train_split`. Again one `test_train_split` is present in the original code and there would have been a,

```
import pandas as pd
from sklearn.model_selection import train_test_split

X = # a dataframe containing the features remove the target y (below)
y = # a same dataframe but a single column which is the target
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=1)

```

The above will give a training proportion at 80% and test proportion at 20%. It's called an 80:20 split. That code will have been in the original code for the first calculation, I know that because the results are pretty standard for a routine ML `test_train_split` operation. Whats happening is the training set is being held absolutely separately to the test set. There is 100% compartmentalisation between the two data sets. Thats critical because if they got confused the calculation falls over, you will be testing the exact data thats just been trained. That confusion would be term "data leakage".

In the code you presented what has happened is the same data used for parameterisation is exactly the same data for training (`X_train`, `y_train`). Again it's called data leakage and parameter grids are very sensitive to this.

To correct this when a parameter grid is used it must have its own data set. That is instead of having a `X_train,X_test,y_train,y_test` the `test_train_split` is modified to generate `X_train,X_test,y_train,y_test, X_val, y_val`. The code to generate this is,

```
import pandas as pd
from sklearn.model_selection import train_test_split
 X_train, X_test, y_train, y_test
    = train_test_split(X, y, test_size=0.2, random_state=1)

 X_train, X_val, y_train, y_val
    = train_test_split(X_train, y_train, test_size=0.25, random_state=1)

```

This approach is singly leveraging `test_train_split` by splitting the training data set twice. This example creates a 60:20:20 split. Some say that 70:20:10 split is used for a parameter split (10 for the parameter search). Matter of opinion (you just adjust the proportions).

**Thus**

1. the grid search is strictly for the parameter split (X\_val, y\_val)
2. and those parameters at then using in training split (X\_train, y\_train)
3. the testing is then performed (X\_test, y\_test)

What this process does is ensure no part of `X_val`, `y_val` is ever seen in the training set. If that happens again - data leakage. Here only the optimised parameters are seen by the training set.

---

**Summary** Introducing validation parameter split into the calculation isn't trivial and there is nothing wrong with your original result. Thus you don't necessarily need it. Generally speaking 80% accuracy is fine, yours is 1% lower (can be overlooked).

Parameter tuning will increase the accuracy score but it will be 85% rather than 79% and if it's really gone well then it might be 90%. What I'm saying is it will improve things, but will not absolutely solve the accuracy issue, it's just a general and accepted issue within ML.

---

Just to finally add there are biological reasons in cancer prediction why accuracy is a key parameter, but thats a different question.

---

Comments: causes lower scores

What I'm saying the correct approach. I recognise that **Geek for Geeks** uses your approach - but they are wrong.

A separate validation split will certainly give result in lower scores than `grid_search.fit(X_train,y_train)` thats inevitable, but it is correct in stopping overfitting.

If the validation split causes its lower accuracy than `0.8` I'd look at `gammas = ... 0.0001,0.001,0.1 ...` because you have missed `0.01` and that will affect your model. It shouldn't give a lower score than the default parameters - it's not possible and something is wrong with the code and application of parameters.

Once you've the initial parameters you can refine the parameter search around the initial values - nothing wrong with that and that can even be coded (automated). The only minor issue is reducing the size of the training split (60% is a bit low). That can cause issues in which case the split is 70:20:10, i.e. 10 for validation.

---
Tags: machine-learning, modelling, microbial-genomics, microbiome
---


# Question
Title: Interpreting GFA graph visualized in Bandage
I assembled Nanopore sequenced reads with Flye and visualized the GFA graph in Bandage but I don't really know how to interpret the result. For context, this is a yeast (DNA) genome. My goal is to figure out areas of repeats, ploidy level, and overall structure of the strain. So far I understand that the bubbles are repetitive areas but I mainly want to know how to interpret a graph like this.

I'm new to interpreting these graphs and I couldn't find any sort of guide online so if someone has any resources they'd be willing to share or if they could help interpret this graph I would really appreciate it!

# Answer
I would suggest making a dot plot of your assembly against a reference (related) yeast genome assembly.

I'd suggest also running QUAST against the yeast reference to get a finer grained idea of the comparison.

Bandage is more of an assembly QC tool for understanding whether there are issues with the assembly. It will not tell you anything about ploidy or that kind of thing.

Off the top of my head, I find it worrisome that 90% of your sequence is in a single scaffold/contig. Presumably there's more than a single chromosome. It does look like the repeats are messing you up.

I'd suggest using hifiasm instead and seeing if it makes different decisions.
> 1 votes

---
Tags: assembly
---


# Question
Title: scRNA: What are good dimensionality reduction/clustering parameters to get biologically plausible groupings?
I've got a moderately large set of PBMCs, over 1M cells. That means I can't easily do a grid search of dimensionality reduction/clustering parameters/methods. Some examples results I'm getting with scanpy default parameters include PCA -- a big, dense grouping with lobes radiating outwards, and UMAP -- "better" groupings in some sense, but many many tiny clusters of cells, more than are biologically plausible (hundreds) as well.

Basically I'm looking for best dimred/clustering practices for such a dataset, other than what's on offer from scanpy/Seurat. My metrics such as they are are that the major cell types should be clearly separated, and that the overall plot have some biological plausibility in terms of number of visually distinct clusters.

I'm aware of the Xiang et al Frontier Genetics 2021 paper, but the number of cells their example datasets use is much smaller AFAIK. They concluded that tSNE is the best choice of method according to their metrics, though UMAP also has some stability advantages. The Chari and Pachter PLOS Comp. Bio. 2023 paper raises some good points (a highly dimensional space can't really accurately have its relations preserved in 2 dimensions, and there are other issues of noise, interpretability, and so on. They recommend clustering from the higher dimensional data, which I do agree with. A different paper rec could be great.

PS: SCVI seems like a good option, I'm playing with that at the moment.

# Answer
> 0 votes
**Aim** It's important to define what all these calculations are doing: they're attempting to extract maximum "multi-dimensional" variance from a data set. The key thing here is that the methodology is old, but it is very well established (long before RNAseq ever existed). I get the point empirical stats is boring, but it's essential to consider whats the underbelly of the methods being used.

To answer the question directly

1. PCA+tSNE is worth considering
2. SCVI-tools is in a different league, and totalIV looks very interesting (what I'm tentatively recommending).

*Note, the commercial sector is often more advanced than public sector: I developed ML methods which intersect it with classical approaches for this basic style of calculation. The scvi-tools however is "full-on" via a DL approach.*

*Note2, I've not read the papers cited except the scvi-tools documentation*

**Data set size** I don't do PMBC calculations - I do perform PMBC *post-hoc* analytics BTW. I personally don't see the data set size as being prohibitive for an unsupervised learning calculation. Pyscan should be able to handle this.

**PCA** is the goto approach stage 1 - but there's a caveat, which is what you're hinting at. First, it's how much variance is trapped within each component. Scanpy will permit this calculation. What you want is \>70% variance, classical stats for example PC1+PC2+PC3+PC4 \> 70%. However, from past posts on this site then for PMBCs this might be unattainable, like no-can- exceed the variance threshold for lots and lots of PCs. That is an issue.

To go through your papers

> Chari and Pachter PLOS Comp. Bio. 2023

Sure, I've not read the paper but on your report, simply taking PC1 and PC2 on bivariate axes is not cool at all. Very little of the variance can be represented. Nobody disputes this. Even if 10 PCs trapped \> 70% of the variance - thats a lot of bivariate plots to assess, e.g. PC1 vs PC3, PC1 vs PC4 etc ... At least you've trapped a good proportion of the variance.

The problem becomes what happens when very little of the variance is captured. This is where Xiang *et al* Frontier Genetics 2021 are likely to be coming from. It's better at that point to use tSNE, but its not ideal.

**Whats the solution?**

**Approach 1** In classical stats (not in RNAseq stuff) the output of PCA would be placed into tSNE and this will cause the variance spread across large numbers of PCs to be trapped by tSNE. This will give you nice clusters. I've not seen this approach at all in RNAseq - it is a common approach. The assumption is that loads and loads of PCs from the PCA trap sufficient amounts of the variance and tSNE will nicely represent that in 2 dimensions. It's possible that is what Xiang *et al* Frontier Genetics 2021 have done, but they may have simply gone straight to tSNE (see note2).

**Approach 2** scvi-tool is impressive and I would recommend totalVI for dimensionality reduction, which they place straight into scanpy.

The method and usage (not the robustness) for totalIV scvi-tools is described here. They are deploying a "hidden Markov"/latent variable style model. See the note below on pss. They've used "CITE-seq RNA", which I don't exactly know what that is but hopefully that will remain relevant.

* It's important to note scvi-tools incorporates lots of other approaches including semi-supervised learning (is that scANVI???), transfer-learning (assume its CNN-style calculations), a standard ANN type model (I think thats scVI). I'm simply proposing what I would guess would work in your case.
* It is worth noting there are imperfections with the authors' methods - in general. Computer vision is a useful example, this used very similar approaches and it took years before serious levels of accuracy were achieved. Trapping absolute accuracy within a single deep learning model "off-the-cuff" I suggest is unlikely.

Semi-supervised learning its basically unsupervised learning but quantifying the weights of each feature - you can use it but I have never got good results from semi-supervised learning methods in general, this might be different though.

---

**ps.** If I've time (i.e. might not happen) I'll pop the link into the related question, where the low amounts of variance between PCs, and how to optimise that is partially answered.

**pss.** Overall scvi-kit for the models outlined need published demonstration they are robust. This part of their work appears partially missing, basically there appears a lot more functionality than published papers. I assume their other stuff is working its way through peer-review. *The totalIV approach: I performed a brief search but was unable to find anything published* \- whilst I'm recommending it - **it appears untested**. What I believe is they've stolen a march on the rest of the field in general, so I'm making a best guess.

---
Tags: rna-seq, single-cell, clustering, pca, dimensionality-reduction
---