# Question
Title: Making a union data from different list of genes
I have two groups of patients

I have a list of driver genes

Clearly each group have mutation in a subset of these gene some of which are common between two groups and some of which are unique to one group

I want a union count matrix of two groups

Let's TP53 only has been mutated in one group, I want TP53 in union matrix putting 0 where absent in one group

This is how my data look

```
> dput(d1)
structure(c(1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), .Dim = c(31L, 26L), .Dimnames = list(
    c("TP53", "TSHZ3", "ARID2", "CDKN2A", "NOTCH1", "APC", "ARID1A",
    "EPHA3", "RNF43", "SCN3A", "SMAD4", "ACVR2A", "ARID1B", "BRAF",
    "CCDC102B", "CCNE1", "CDH1", "CHD4", "EGFR", "EPHA2", "ERBB2",
    "FBXW7", "KDM6A", "LIN7A", "MAP2K7", "MET", "MSH3", "PIK3CA",
    "PTEN", "RASA1", "SLIT2"), c("LP6008336-DNA_G01", "LP6008334-DNA_B02",
    "LP6008334-DNA_A04", "s24", "s15", "LP6008334-DNA_C02", "LP6008460-DNA_D01",
    "LP6005500-DNA_D03", "LP6008337-DNA_A07", "LP6007600", "LP6008337-DNA_H06",
    "LP6005334-DNA_H01", "LP6008460-DNA_A04", "LP6008334-DNA_A03",
    "LP6008334-DNA_D02", "LP6008336-DNA_H01", "s30", "s86", "s59",
    "LP6008336-DNA_F02", "s67", "LP6008460-DNA_F02", "s94", "s18",
    "LP6008460-DNA_G03", "s80")))
> dput(d2)
structure(c(1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L,
2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 2L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L,
0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 2L, 2L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L), .Dim = c(33L, 35L), .Dimnames = list(
    c("TP53", "NAV3", "CDKN2A", "SMAD4", "KCNQ3", "ARID1A", "ARID2",
    "LRRK2", "EPHA3", "KRAS", "LIN7A", "PIK3CA", "PTEN", "RNF43",
    "NOTCH1", "PBRM1", "SCN3A", "SLIT2", "TSHZ3", "APC", "ARID1B",
    "CCDC102B", "CTNNB1", "EGFR", "ERBB2", "GPATCH8", "MSH3",
    "NIPBL", "POLQ", "SIN3A", "SMARCA4", "STK11", "ZFHX3"), c("LP6005500-DNA_H03",
    "LP6008031-DNA_B02", "LP6008269-DNA_E10", "LP6005334-DNA_B02",
    "LP6008269-DNA_H05", "LP6008031-DNA_D02", "LP6005500-DNA_G01",
    "LP6005334-DNA_H03", "LP2000329-DNA_A01", "LP6005690-DNA_F02",
    "LP6008337-DNA_B07", "LP6008031-DNA_B03", "LP6005500-DNA_C03",
    "LP6005690-DNA_A02", "LP6005690-DNA_E01", "LP6008334-DNA_B01",
    "LP6008051-DNA_B01", "LP6005409-DNA_B02", "LP6005409-DNA_F01",
    "LP6005935-DNA_D02", "LP6005935-DNA_F01", "LP6008221-DNA_C03",
    "LP6008280-DNA_D02", "LP6008337-DNA_F06", "LP6008202-DNA_F02",
    "LP6008280-DNA_B02", "LP6005334-DNA_C02", "LP6008031-DNA_E01",
    "SM-4AX87", "LP6005409-DNA_E02", "LP6008031-DNA_D01", "LP6008202-DNA_A03",
    "LP6008051-DNA_A01", "LP2000325-DNA_A01", "LP6005409-DNA_H01"
    )))
>

```

I have tried

```
plyr::rbind.fill.matrix(d1,d2)

rbind(d1,d2,all=TRUE)

```

But the outputs seem crab while need a trustful data for machine learning

# Answer
Maybe something like this:

```
df = merge(
data.frame(gene=rownames(d1),d1),
data.frame(gene=rownames(d2),d2),all=TRUE)

```

If you want to replace NA with 0:

```
df[is.na(df)] = 0

```
> 1 votes

---
Tags: r, merge
---


# Question
Title: Find rsIDs for GRCh37 SNPs and rsIDs for GRCh38 SNPs and compare for overlap
I have a long list of variant IDs that I generated as a result of research done one GRCh37 genomes (e.g. `13_28025615_G_C_b37`). I want to get their rsIDs and compare them to a newly-obtained list of GRCh38 rsIDs that have also been retrieved from similarly formatted variant IDs (e.g. `chr22_50473581_G_A_b38`). How would I go about doing this? I've tried installing `biomaRt` as a package but I can't quite make heads or tails of it and it's not clear to me that it would support both versions of genome builds.

# Answer
> 3 votes
> it's not clear to me that it would support both versions of genome builds.

You can specify what version of the database you are querying.

I'd start with biomart on the ensembl website, see if you can get it to do what you want, then work on getting biomaRt to do what you want.

# Answer
> 2 votes
rsIDs change with the dbSNP builds, not with the reference genome builds. rsIDs do not depend on a reference genome. They point to a specified locus regardless of the differences in genomic assemblies. rsIDs numbersare assigned by dbSNP database maintainers. The numbers may be deleted from database (e.g. when a submitter withdraw data) and may change on merges (e.g. when multiple rsIDs were assigned at the same genomic location).

Specifying the location by rsID is stable, and can be only affected by an eventual merge or deletion, if any, with a new dbSNP database build. It cannot be affected by a reference genome build.

However, a base position (offset) of a locus within a chromosome may change between the builds of the reference genome. So, when specifying the location by the chromosome and base position, the reference genome build should be also given, to locate an appropriate rsID.

You can download dbSNP database files (in VCF format) of the database build that was actual at a time when a specific build of a reference genome was also actual, but these files are huge. A possible place to download is https://ftp.ncbi.nih.gov/snp/pre\_build152/organisms/ or https://ftp.ncbi.nih.gov/snp/archive/

However, these files are arranged by the dbSNP database build (e.g. 152, 153, 154), not by the reference genome build (e.g. 37 or 38).

If your list is relatively small, you can fetch the data rs numbers via dbSNP API by chromosome and location. The API is available at https://www.ncbi.nlm.nih.gov/snp/docs/eutils\_help/

For instance, do “advanced search” of dbSNP at https://www.ncbi.nlm.nih.gov/snp/advanced When it comes to search by the base position, you can only specify positions for GRCh37 (called “Base Position Previous”) and GRCh38 (“Base Position”). This URL will help you construct the request string that you will further use for the dbSNP API requests.

You may also need to add “NOT MergedRS\[All Fields\]” to exclude rsIDs that have already been merged. I wish I knew what to type instead of “All Fields” here to refine the search, but the documentation is very sparse.

In your examples, for `13_28025615_G_C_b37` specify Chromosome=13 and Base Position Previous=28025615 in the “Advanced Search”.

As about the API, for `13_28025615_G_C_b37` query https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=13\[Chromosome\]+AND+28025615\[Base+Position+Previous\]+NOT+MergedRS\[All+Fields\]

and for `chr22_50473581_G_A_b38` query https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=22\[Chromosome\]+AND+50473581\[Base+Position\]+NOT+MergedRS\[All+Fields\]

The first rsID that you will get will be the main current rsID, the others (if any) will be the merged ones.

Here is a small perl example on how to retrieve data from the dbSNP. The fetch\_GRCh37 function returns the rsID for a given chromosome and base position. For GRCh38, remove the “+Previous” from the URL. You can call the fetch\_GRCh37 in a script multiple times. It’s better that you get rsIDs for your both GRCh37 and GRCh38 data to make sure that you refer to the same rsIDs that are actual after eventual merges.

Please make sure to read and comply to the requirements before using the API.

A sample Perl script snippet follows:

```
use XML::LibXML;
use LWP::UserAgent;

sub fetch_GRCh37 {
        my $chromosome = $_[0];
        my $basepos = $_[1];
        my $url = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term='.$chromosome.'[Chromosome]+AND+'.$basepos.'[Base+Position+Previous]+NOT+MergedRS[All+Fields]';
        my $UserAgent = new LWP::UserAgent ();
        my $GET = $UserAgent->get($url)->content;
        my $dom = XML::LibXML->load_xml(string => $GET);
        my $numberofnodes = int($dom->findnodes('eSearchResult/Count')->to_literal());
        if ($numberofnodes > 0) {foreach my $title ($dom->findnodes('eSearchResult/IdList/Id')){return("rs".$title->to_literal());}}
}

print fetch_GRCh37(13, 28025615);

```

---
Tags: variation, biomart
---


# Question
Title: Convert rs ID of one hg build to rs IDs of another build
I have a list of dbSNP rsIDs for GRCh37 and I want to convert them to the equivalent IDs in GRCh38. This is using the most recent dbSNP build (150 as of the time of this post). Is there any ID mapping available? If not, are there any tools I can use?

# Answer
> 2 votes
Did you try the liftOver tool of UCSC? You'll need a BED file with your SNPs coordinates for this.

You can also do this in `R` with `rtracklayer`.

```
library(rtracklayer)
?liftOver

```

# Answer
> 5 votes
You can assume that the overwhelming majority of rsIDs are the same between GRCh37 and GRCh38 (they're semi-stable IDs). There are, however, a number of rsIDs that are present only in GRCh37, which you can find here. Note that the format of this file is a bit strange, it's `chromosome|position|ID|weight`, where `position` is sometimes empty, `weight` is typically 1, and `ID` is the rsID without the `rs` prefix.

Also, some rsIDs have a different strand. This is less likely to be an issue for you, but it's there in case you need it.

# Answer
> 1 votes
See also CrossMap: http://crossmap.sourceforge.net/

20, 21, 22, 23, 24, 30

# Answer
> 0 votes
These are only the base positions that have changed between GRCh37 and GRCh38, not the IDs. For example if a study (e.g. found via pubmed) refers to a particular rsID, the genome build does not matter.

rsIDs do not depend on a reference genome. They point to a specified locus regardless of the differences in genomic assemblies. As a result, you cannot convert rsIDs of one genome build to rsIDs of another build. And there is no such mapping tool available, even if they pretend to be. rsIDs are assigned by dbSNP maintainers, and are stable. The numbers may only change on merge, where multiple numbers were assigned at the same genomic location, and the higher number gets merged into a lower number by the dbSNP maintainers, unless there is a very well-published, notable rsID mentioned in the press, in this case that merged number is considered the mentioned one, not the smaller one.

> RefSNP number is the stable accession regardless of the differences in genomic assemblies.

For more information go to the “RefSNP Number Stability” section at https://www.ncbi.nlm.nih.gov/books/NBK21088/

The rsIDs may also be deleted in some cases see the following links on the explanation:

Thus, specifying the location by rsID is stable, and can be only affected by an eventual merge or delete from the dbSNP database, if any, not by a new genome build. You can only “update” the rsIDs from time to time to make sure they all are the main rsID, not those larger numbers already merged to a smaller one. It should be noted that the base positions (offsets) of the same genomic location may change between the builds of the reference genome. And they did surely change between GRCh37 and GRCh38.

If your list is relatively small, you can get most recent rsIDs for your chromosome and offset via dbSNP API. It is available at https://www.ncbi.nlm.nih.gov/snp/docs/eutils\_help/

The query for the API is the same as for the “advanced search” of dbSNP at https://www.ncbi.nlm.nih.gov/snp/advanced

If you specify data in a “Base Position” line, it will look for GRCh38 base position. If you need to specify the base position for a GRCh37, use “Base Position Previous” line. You may also need to add `NOT MergedRS[All Fields]` to exclude rsIDs that have already been merged.

For example, to get the rsID for a chromosome 1 GRCh37 base position 249222527, the API query URL is the following https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=1\[Chromosome\]+AND+249222527\[Base+Position+Previous\]+NOT+MergedRS\[All+Fields\]

You can do queries from a script if the number of requests is relatively small. You may find a perl script, for example, that does the job.

Even though the rsID numbers can be merged or deleted, and the changes are reflected in a new dbSNP database build, given the fact that the changes in numbers do not depend on the genome build, is better to always rely to the rsID number from the latest dbSNP database build.

The changes in rsID numbers are changed due to merges or deletions between dbSNP builds relatively infrequently, given the total number of the IDs in the database.

---
Tags: snp, plink, phylogenetics
---


title: GWAS MAC filter Interpretation

text:
# Question
Title: GWAS MAC filter Interpretation
I am performing a GWAS analysis and try to understand the influence of the minor allele count filter.Setting the filter to 1 % gave me this plot and I am confused about the same -log10 pvalues around ~ 3.8 which shows a strange pattern (many dots in one line with same value).

If I increase the filter to 6 % the dots in a line disappear but I loose also some significant SNPs:

Sample are about 77000 SNP positions (barley), unimputed from GWAS (all NaNs removed). GWAS was performed with 200 phenotype samples.

How to interpret the dots in a line?

# Answer
> 1 votes
This looks like genotyping failed for some samples. It's unclear from the op whether the genotyping data was controlled for quality, but I would recommend performing sample QC.

On a side note: in human gwas with imputed snps, I would expect to see clusters of significant snps as a sign that the significant snps are not an artifact. In the op there are some single significant snps. This might be a sign that those significant snps are artifacts. There are also some snps that form clusters and therefore look more promising.

---
Tags: python, gwas
---


# Question
Title: the variation between treatments is less than the variation between replicates in RNA-seq data
I have a set of `RNA-seq` samples from targeting different proteins in a complex with `siRNAs`. However, the `PCA` plot shows that **the distance between the two replicates of each treatment is more than the distance between the treatments**. So, there is **no clear border between treatments** to cluster them separately.

If I go through my usual RNA-seq data analysis pipeline using `EdgeR` normalization and `Limma/Voom` DEG analysis, I get a lot of DEGs. However, when I try to visualize them using a heatmap, the replicates do not cluster together, as expected.

What should I do in this situation? I mean should I trust the DEGs, or try to still extract co-expressed gene modules from this data?

my code for normalization:

```
myDGEList <- DGEList(Txi_gene$counts)
log2.cpm <- cpm(myDGEList, log=TRUE)

cpm <- cpm(myDGEList)

# prepare your grouping matrix
group <- factor(targets$group)
design <- model.matrix(~group)

# filtering
keepers <- filterByExpr(cpm, design)
myDGEList.filtered <- myDGEList[keepers,]
log2.cpm.filtered <- cpm(myDGEList.filtered, log=TRUE)

myDGEList.filtered.norm <- calcNormFactors(myDGEList.filtered, method = "TMM")
log2.cpm.filtered.norm <- cpm(myDGEList.filtered.norm, log=TRUE)

```

my code for PCA plot:

```
pca.res <- prcomp(t(log2.cpm.filtered.norm), scale.=F, retx=T)
pc.var <- pca.res$sdev^2
pc.per <- round(pc.var/sum(pc.var)*100, 1)
pca.res.df <- as_tibble(pca.res$x)

# PC1 vs PC2
pca.plot <- ggplot(pca.res.df) +
  aes(x=PC1, y=PC2, label=sampleLabels, color = group, shape = batch) +
  geom_point(size=4) +
  xlab(paste0("PC1 (",pc.per[1],"%",")")) +
  ylab(paste0("PC2 (",pc.per[2],"%",")")) +
  labs(title="PCA plot",
       caption=paste0("produced on ", Sys.time())) +
  coord_fixed() +
  theme_bw()

ggplotly(pca.plot)

```

the PCA plot figure. The legend is cropped out for privacy reasons. Colors represent replicates. Shapes represent batch effect which was has been fixed using `combat_seq` on raw-data.

# Answer
> 1 votes
If you got DEGs with statistically significant p-values after multiple testing following a valid pipeline, **I'd believe that**, even if looking at the first few PCs doesn't look promising.

---
Tags: rna-seq, deseq2, edger, pca, limma
---


# Question
Title: Analyzing transcripts for non-coding RNA
From a RNA seq experiment we look for the counts of transcripts of tRNA but we are also interested of other types of RNA we can find. Because of the protocol used for the library prep we can't analyse miRNA but we can analyse ncRNA. I never analyzed transcripts for ncRNA and I have some questions:

* The name of the transcripts are like ENST00000611544.1, ENST00000387069.1, since they are non coding, are they linked to a gene? Should I rename the name of each ncRNA with a gene name or working with these long names is ok?
* When analyzing transcripts of RNA I did enrichment analysis to see what pathways are more active. But when analyzing transcripts for ncRNA is possible to do something similar to enrichment analysis? How can analyse the functions of the ncRNA that we get from the study?

Thank you very much in advance

# Answer
> 2 votes
ncRNAs are by itself not special. Sure they might be more unstable and less abundant than other transcripts but for the analysis it is more or less the same. You should take your transcript level counts and aggregate them to the gene level. The tximport package comes to mind. Every transcript (given it is annotated) is linked to a gene. ENST00000611544 is connected to the U2 gene but these snRNAs are also short and I doubt you see them in normal RNA-seq.

In the end you can process data as any RNA-seq. Quantify your data against a transcriptome, e.g. with `salmon`, then aggregate transcript level abundance estimates (so basically counts) to the gene level and see whether your genes of interest have counts you can work with or whether you mostly have zeros. A transcriptome can be obtaine e.g. from GENCODE. The available annotation files will tell you which genes these transcripts originate from. Check for the GTF files.

Towards enrichment analysis, well it all comes down to the annotations you have. If the ncRNAs are annotated in any pathway then you can perform standard analysis with tools such as goseq or gprofiler2. If they are not annotated, well then there is no information available.

---
Tags: rna-seq
---


# Question
Title: Where can I find the genotype frequencies for the SNP's that are tested at 23andMe, Ancestry and FTDNA/MyHeritage?
I want to add the frequencies of the genotypes to a baseline of raw DNA data I've got from the DNA testing companies 23andMe, Ancestry and MyHeritage (tested at all). All in build 37 human genome.

Where can I find that information? Checked HapMap project to no avail, can't align with the rsId or positions used in the text files.

Thanks in advance for your help

# Answer
gnomad is probably the most comprehensive allele frequency database that exists at the moment.

Depending on how many variants you need to query, you can either manually input them into the search bar, or if you have a lot more, then you can also directly download the vcfs under the downloads tab.

The default for V2 gnomad set is on the gr37 build.
> 2 votes

# Answer
The largest number of samples is at the dbSNP database. To view frequencies of a single rsID in a web browser, for example, rs11950646, use the following URL:

https://www.ncbi.nlm.nih.gov/snp/rs11950646#frequency\_tab

The dbSNP has more samples to calculate the frequency than gnomad, at least for this rsID.

If you need to automatically process multiple rsIDs, as you said, for the whole genotype file that you got from the companies like 23andMe, you can download a VCF with allele frequency data from https://ftp.ncbi.nih.gov/snp/population\_frequency/latest\_release/ and write a script to automatically process the data and make the output in a format you want. If you only need frequency of some rsIDs within a single chromosome, you can download the whole database for that chromosome from https://ftp.ncbi.nih.gov/snp/latest\_release/JSON/ and will find the data in the JSON files, each file for each chromosome (plus the mtDNA). Look for the "frequency" field.

If you are not engaged to write scripts, and just wish to interactively browse various rsIDs in context of a frequency, you can upload your genotype file to selfdecode.com - it is a convenient, albeit a paid service.
> 2 votes

---
Tags: snp, public-databases
---


title: PDB residue numbering

text:
# Question
Title: PDB residue numbering
The residue numbering in PDB (in my understanding) coincides with UniProtKB align. However, this might not be consistent with that in PDB fasta file. To complicate it further, UniProtKB align can be a few mutations away from the PDB sequence. Disordered/unidentified residues in structures also create gaps in numbering.

Is there a way/suggestion to create a matched residue numbering between sequence and structure? For example

```
12345678 (residue numbering)
ATPWQMSQ (sequence)
---W-MSQ (PDB residues) "-" for missing residue in structure

```

# Answer
> 2 votes
Unfortunately, PDB numbering can be one of different schemes. The main three:

* Uniprot numbering (canonical isoform)
* Construct numbering with or without tag removal
* Start from 1 for what has density

However, the are many many corner cases. In Proteopedia there is a lovely article about weird PDB numbering.

The best way to find the details about what part of the sequence is covered is to use the Swissmodel API. It's actually better than either the PDBe and RCSB PDB APIs. Swissmodel is a homology modelling service that maintains close threaded models (\>30% ID), but also has PDB data. Here is a Python3 method I use, (`provider='pdb'`) is what is required here).

```
def get_data(uniprot: str, provider='swissmodel') -> List[dict]:
        """
        Gets data from Swissmodel. Convert each entry with convert_SM_entry
        """
    assert provider in ('swissmodel', 'pdb'), 'Provider has to be swissmodel or pdb'
        url = f'https://swissmodel.expasy.org/repository/uniprot/{uniprot.strip()}.json?provider={provider.strip()}'
        data = requests.get(url).json()
        return data['result']['structures']  # [0] # from to coordinates description

```

The returned dictionary contains a start and end of the uniprot \<--\> PDB. It does not however, say where the missing density residues are. But in *deposited* structures the PDB numbering skips where these residues ought to be... So might not be an important problem.

Another database of interest is SIFTS.

To make a multiple sequence alignment of a non-deposited PDB, you need to extract the sequence and do the alignment here is a snippet using pymol module (from conda):

```
import pymol2

pdbfile = '' #your pdb file name..
selector = 'chain A' # or whatever chain you want
with pymol2.PyMOL() as pymol:
    pymol.cmd.load(pdbfile)
    fasta = pymol.cmd.get_fastastr(selector)

seq = ''.join(fasta.replace('?','').split('\n')[1:]) # remove header and non AA residues

```

Biopython has a dozen sequence alignment tools. But the original one still works fine

```
from Bio import pairwise2
from Bio.Seq import Seq

alignment = pairwise2.align.globalxx(''.join(seq).replace('?',''), "ACQQG")[0]

```

`alignments` is a tuple of the MSA sequences.

If you want to do the latter without coding. Extract the sequence from the PDB using PyMOL the normal application and load the pdb and then `print cmd.get_fastastr('chain A')` in the command line part of the GUI. Copy this. Go to Muscle or other onliner aligner and align!

---
Tags: proteins
---


# Question
Title: RNAseq biological replicates not clustering in PCA plots
I have RNAseq data from 4 samples with 3 biological replicates per sample. I am currently trying to do the differential expression analysis with DESeq2 but the biological replicates will not cluster together when I make the PCA plot or correlation heatmap. This is my first time with RNASeq analysis and so am not sure what the best route forward is? I would like to avoid repeating the experiment with new samples if possible!

My pipeline prior to DESeq2 was the following:

FastQC quality check -\> Trimmomatic -\> Kallisto

I used tximport to convert kallisto files into suitable format for DESeq2

# Answer
> 1 votes
PC1 is 81% of the variance?

This PCA plot confirms that different tissues are different. You probably already knew that. I'd make more PCA plots that are tissue specific. That will be more informative than these.

Personally, I'd also not do DESeq on all of these samples together, unless the goal of your experiment really is to learn the differences between these two tissues.

# Answer
> 1 votes
You used Kallisto for alignment. I think Kallisto reports TPM values, are you using this value? DEseq2 uses count data, so I am not sure whether these two methods are compatible.

Also, I agree with previous answers that your PCA actually looks OK. One possible way to improve is to choose top variable genes. For example, you can try top 3,000, 5,000, 7,000 genes and so on. The idea is that for the genes that do not show much variation between samples, including them in PCA may just introduce noise. You can also try to color samples in your PCA by some other variables, like batch, sequencing depth and so on to trouble-shooting.

# Answer
> 0 votes
Clustering of replicates looks decent enough to me, so you should be abl to push ahead, but I agree the tissues are grouping, which could mask any differences based on sex or genotype.

You might consider the EdgeR package for DE analysis here. It allows for flexibility when making complex comparisons while accounting for tissue/batch effects. I've had good luck with using it to compare across batch effects from complex experiments like this.

---
Tags: rna-seq, deseq2, clustering, pca
---


# Question
Title: plink: --update-name vs. editing the BIM
I am working with data from the HGDP public files. These are massive VCFs from whole genome sequencing. To begin with, I converted these to BEDs using plink (1.9). The next problem is that these files use rsid names instead of chrpos names. I would like to change to chrpos. I am wondering whether one can simply edit the BIMs or whether I have to do use plink's `--update-name` call. The problem with the latter is that it doesn't support duplicates. The HGDP files are full of duplicate variant names because many of the variants don't have rsid's and are just named ".":

```
$> plink --bfile hgdp_chr21 --update-name updateids --make-bed --out hgdp_chr21_chrpos
PLINK v1.90b6.17 64-bit (28 Apr 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to hgdp_chr21_chrpos.log.
Options in effect:
  --bfile hgdp_chr21
  --make-bed
  --out hgdp_chr21_chrpos
  --update-name updateids

64259 MB RAM detected; reserving 32129 MB for main workspace.
1086522 variants loaded from .bim file.
929 people (0 males, 0 females, 929 ambiguous) loaded from .fam.
Ambiguous sex IDs written to hgdp_chr21_chrpos.nosex .
Error: Duplicate ID '.'.

```

Thus, one gets an error when trying `--update-name` and I am wondering if I can simply edit the BIM files to be chrpos instead of rsid. I can't find an answer to this anywhere. It seems that if the BED files don't contain a reference to the variant IDs, but only to the row in the BIM file, then this approach should be possible.

Here's the current BIM file:

```
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38.bim
1   rs1331781659    0   10153   G   A
1   rs1179689498    0   10163   C   T
1   rs201694901 0   10180   C   T
1   rs199706086 0   10250   C   A
1   rs111200574 0   10257   C   A
1   rs145427775 0   10291   T   C
1   rs868413313 0   10297   T   C
1   rs112750067 0   10327   C   T
1   rs1035171912    0   10330   A   C
1   rs868804735 0   10333   T   C

```

As can be seen, very few lines have the "." variant name. In fact, in the first 100 lines there are 0.

# Answer
> 2 votes
I would reccomend using `bcftools` on the original vcf files before you convert them to plink, to fill in missing IDs using the command:

```
bcftools annotate --set-id +'%CHROM\_%POS\_%REF\_%FIRST_ALT' file.vcf

```

This means you won't have any duplicate ID names and plink should be happy about that. It will preserve the positional information in the .bim file as well.

# Answer
> 2 votes
I found this 2019 question on biostars in which user zx8754 mentions that plink2 has a command for this purpose `--set-all-var-ids`, from the plink2 docs:

> Whole-exome and whole-genome sequencing results frequently contain variants which have not been assigned standard IDs. If you don't want to throw out all of that data, you'll usually want to assign them chromosome-and-position-based IDs.
>
> --set-missing-var-ids (which just replaces missing IDs) and --set-all-var-ids (which overwrites everything) provide one way to do this. The parameter taken by these flags is a special template string, with a '@' where the chromosome code should go, and a '#' where the base-pair position belongs. (Exactly one @ and one # must be present.) For example, given a .pvar file starting with
>
> #CHROM POS ID REF ALT chr1 10583 . G A chr1 886817 . T C chr1 886817 . C CATTTT
>
> "--set-missing-var-ids @:#\[b37\]" would name the first variant 'chr1:10583\[b37\]', the second variant 'chr1:886817\[b37\]'... and the third variant also gets the name 'chr1:886817\[b37\]'.

Thus, in my case, I used:

```
user@desktop:/media/luks8tb1/data/genomics/HGDP$ plink2 --bfile hgdp_hg38 --set-all-var-ids @:# --make-bed --out hgdp_hg38_chrpos
PLINK v2.00a3LM AVX2 Intel (27 Jul 2020)       www.cog-genomics.org/plink/2.0/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to hgdp_hg38_chrpos.log.
Options in effect:
  --bfile hgdp_hg38
  --make-bed
  --out hgdp_hg38_chrpos
  --set-all-var-ids @:#

Start time: Tue Aug 25 19:16:10 2020
64259 MiB RAM detected; reserving 32129 MiB for main workspace.
Using up to 8 compute threads.
929 samples (0 females, 0 males, 929 ambiguous; 929 founders) loaded from
hgdp_hg38.fam.
56475980 variants loaded from hgdp_hg38.bim.
Note: No phenotype data present.
Writing hgdp_hg38_chrpos.fam ... done.
Writing hgdp_hg38_chrpos.bim ... done.
Writing hgdp_hg38_chrpos.bed ... done.
End time: Tue Aug 25 19:17:41 2020

```

Results, BIM:

```
emil@test-desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38_chrpos.bim
1   1:10153 0   10153   G   A
1   1:10163 0   10163   C   T
1   1:10180 0   10180   C   T
1   1:10250 0   10250   C   A
1   1:10257 0   10257   C   A
1   1:10291 0   10291   T   C
1   1:10297 0   10297   T   C
1   1:10327 0   10327   C   T
1   1:10330 0   10330   A   C
1   1:10333 0   10333   T   C

```

---
Tags: bed, plink
---


# Question
Title: how to estimate phylogeny without DNA sequences?
I am new to bioinformatics, I am reading ‘Analysis of Phylogenetics Second Edition and Evolution with R’ from Emmanuel Paradis.

I can create phylogeny from DNA sequences, by first calculating the distances and then using (ie nj() function).

My question is: if instead of DNA I have trait data (ie height and weight of species), how can I estimate a phylogeny? (is the same process valid? somehow calculate distances and pass that to nj()?)

# Answer
> 2 votes
Sequence data is basically always better and at this point in history it is de rigueur for any serious study, but originally phylogenies were all based on (ideally discrete) phenotypic characters, so there's no reason why it can't work in principle. As @Pallie mentions in a comment, the clustering methods in e.g. NJ or UPGMA will work on any set of distance measures, and they were in fact originally developed on phenotypic characters rather than sequence data.

I'd recommend taking a look at chapter 24 of Joe Felsenstein's book on phylogenies, which is digitized here (yes the link looks sketchy, sorry). There is a section on p. 425 "Inferring phylogenies and correlations", which discusses how you can in principle use quantitative characters to infer phylogenies. The theory and math for why this is supposed to work is earlier in the chapter.

In many cases people prefer to code the characters as discrete because it makes it easier to work with. But there are also methods that work with the raw quantitative characters.

In R, you should be able to do something like this:

```
# some simulated data, species are rows and characters are columns
> data
      char1 char2 char3
spec1     0     3     6
spec2     1     4     7
spec3     2     5     8

> library(ape)
> tr = nj(dist(data))
> plot(tr)

```

---
Tags: phylogenetics, phylogeny
---


# Question
Title: How to annotate gene length to a list of gene symbols using UCSC data?
I have a list of HGNC gene symbols, I am looking to get the gene length of each gene. Although I also describe these genes with lots of UCSC datasets as features, so I am wondering if there is a dataset in UCSC I can use to also get gene length from?

I've been looking in the data that is downloadable from UCSC table browser (I've been aiming to find start and ends for each gene to subtract to get gene length) but there are a lot of files and I'm not sure which dataset to take from which will also match to my HGNC gene symbols.

# Answer
> 2 votes
You could do a few command-line operations to answer this question. This assumes the use of `hg38` assembly.

First, get a list of genes from GENCODE:

```
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

```

Then use `grep` to filter those genes with your list of HGNC symbols:

```
$ grep -wFf hgnc_symbols.txt genes.bed > filtered_genes.bed

```

You can run this through `awk` to get lengths:

```
$ awk -vFS="\t" '{ print $3-$2 }' filtered_genes.bed > filtered_gene_lengths.txt

```

---
Tags: gene, genomics, ucsc
---


# Question
Title: probability of an entire sequence, not just one value
I am reading "`Statistical methods in bioinformatics`"

I came across this formula:

As explained in the book **O** are the observations:

My question is why the author only specified $prob(O)$ for a sequence of length=1 ?, I assume he is only talking about one observation, but generally we have a sequence $O\_1 ,O\_2,O\_3,...$.

I assume that $Prob(sequence)= P(O\_1).P(O\_2)...$ but I am not sure.

# Answer
The author writes in the quote that $O$ is to be taken for the whole set of positions $O\_{i}$.

In the extracts that you provide there is no reason to believe that it is only for a sequence of length 1. The author writes "the observed sequence $O = O\_1, O\_2...$"

Thus as you say $Pr(sequence) = Pr(O) = Pr(O\_1) * Pr(O\_2)...$ according to my understanding of the quotes.
> 2 votes

---
Tags: ngs, hidden-markov-models
---


# Question
Title: How does DeepVariant construct RGB images from DNA sequences?
DeepVariant is a pipeline to call genetic variants from DNA sequencing data.

A major step, before feeding the CNN, is to translate these DNA sequences into images. It's unclear why and how Google constructs the RGB images from the DNA data. Obviously, DNA is a string over an alphabet with the characters: {A, T, C, G}.

It is even hard to understand how the mapping works based on the source code of their unit tests.

In their figure from the paper: A is Red, C is Green, G is blue, and T is Yellow (G+R), but this is still unclear how they construct the `3xNxN` image.

---

EDIT from google's blog:

> In this article we will show the six channels in a row, but in DeepVariant they are encoded as six layers in the third dimension, giving each tensor a shape of (100, 221, 6) corresponding to (height, width, channels). The variant in question is always in the center of each pileup image, here marked with a small line at the top.
>
> Channels are shown in greyscale below in the following order:
>
> Read base: different intensities represent A, C, G, and T.
>
> Base quality: set by the sequencing machine. White is higher quality.
>
> Mapping quality: set by the aligner. White is higher quality.
>
> Strand of alignment: Black is forward; white is reverse.
>
> Read supports variant: White means the read supports the given alternate allele, grey means it does not.
>
> Base differs from ref: White means the base is different from the reference, dark grey means the base matches the reference.

# Answer
Actually the paper has made it clear how they did it. You just have to read the supplementary materials closer.

> In their figure from the paper: A is Red, C is Green, G is blue, and T is Yellow (G+R), but this is still unclear how they the 3xNxN image.

In RGB, each dimension is an NxN image. Since you have three dimensions, so it's 3xNxN. The red dimension was used to encode the nucleotide bases. The green dimension was used to encode quality scores. Finally, the blue dimension was used to encode the strand information.

> Obviously, DNA is a string over an alphabet with the characters: {A, T, C, G}.

This is easy. Lot's of ways. You could do a one-hot encoding or what DeepVariant used:

```
def get_base_color(base):
    base_to_color = {'A': 250, 'G': 180, 'T': 100, 'C': 30}
    return base_to_color.get(base, 0)

```
> 1 votes

# Answer
Based on the supplementary material as mentioned in a comment by Devon Ryan:

> The second phase of DeepVariant encodes the reference and read support for each candidate variant into an RGB image. The pseudo-code for this component is shown below; it contains all of the key operations to build the image, leaving out for clarity error handling, code to deal with edge cases like when variants occur close to the start or end of the chromosome, and the implementation of non-essential and/or obvious functions.

Here is the main function

```
WIDTH = 221
HEIGHT = 100;
def create_pileup_images(candidate_variants):
 for candidate in candidate_variants:
 for biallelic_variant in split_into_biallelics(candidate):
 start = biallelic_variant.start - (WIDTH-1) / 2
 end = WIDTH - span_start
 ref_bases = reference.get_bases(start, end)
 image = Image(WIDTH, HEIGHT)
 row_i = fill_reference_pixels(ref, image)
 for read in reads.get_overlapping(start, end):
 if row_i < HEIGHT and is_usable_read(read):
 add_read(image, read, row_i)
 row_i += 1
 yield image

def fill_reference_pixels(ref, image):
 for row in range(5):
 for col in range(WIDTH):
 alpha = 0.4
 ref_base = ref[col]
 red = get_base_color(ref_base)
 green = get_quality_color(60) # The reference is high quality
 blue = get_strand_color(True) # The reference is on the positive strand
 image[row, col] = make_pixel(red, green, blue, alpha)
 return 5

```
> 0 votes

---
Tags: python, variant-calling, machine-learning, deep-learning
---


# Question
Title: How to know which ortholog database is the best for my work?
I am working on four data sets of yeast genes and one data set of mouse gene. I need the list of human orthologs of those genes. There are a lot of databases available online which is in use for this purpose

https://questfororthologs.org/orthology\_databases

I don't know which one to choose and why. Can anyone tell me on what basis should I choose a database over any other database? Thank you.

# Answer
> 0 votes
**To recover orthologous sequence for model species**

You should be able to find your orthologs using BioMart and your gene IDs.

There are numerous way to do it, directly on the ensembl website or through API. You can start from these links (quite old though) :

**Which database should I use**

It is broad and depends on the nature of your work. In your case, working with "major model species", **DIOPT** (the first line of the table in Quest For Orthology) has the very feature you need on its front page for going from mouse to human. More eukaryotic/viral oriented database will be used for taxonomic identification, annotation inference, etc.

I hope it helps.

---
Tags: orthologues
---


# Question
Title: Populations genetics and dynamics of bacteria on a Graph
**Disclaimer**:

I'm a physicist and I'm fairly new to Bioinformatics therefore somethings below may not make sense.

---

**My purpose**:

I would like to simulate genetic evolution of bacterial populations, implementing interactions among different bacterial strains on a graph.

My starting point is to generate a graph with a certain topology.

Each vertex $V\_i$ should contain a population of $N\_i$ individuals for which I choose a certain allelic profile from an MLST database. Population numerosity and allelic profile is assigned by random choice. $$\\\\\$$

At starting time each vertex is populated with a "pure" strain of $N\_i$ clones, identified with a randomly choosen allelic profile.

**Dynamics assumptions**

I will assume that:

* Topology does not change with time
* Vertex dynamic flows as a *Moran process*: At time $t$ , one individual is chosen randomly to reproduce and one individual is chosen to die. The same individual can be chosen to reproduce and then die. Thus, an individual has either zero, one, or two descendants. Zero and two with equal probability $p\_0 = p\_2 = (N\_i -1 )/N\_i^2$ , and one with probability $p\_1 = 1 − 2 p\_2$.
* Migration from one Vertex to one another is possible for nearest neighborhood with some fixed probability depending on choosen topology. I will assume that when migration occurs $I$ new individuals enters into $V\_i$ and $O$ individuals exits from $V\_i$, with $I$ and $O$ depending on origin vertex numerosity. Migrants carry allelic frequencies according to population they came from. When migration occurs a new population is established: allelic profile frequencies changes according to that.
* Migration and *Moran process* evolution does not happen simultaneously; This allows me to guarantee that each population performs at least one evolution.

---

**Questions**:

(1) Is there something fundamental that I'm missing?

(2) Moran model ensures that the population size remains constant (apart from migration effects obviously). Is this appropriate in order to simulate dynamic in stationary growing phase or should I introduce some demographic effects?

(3) Which could be a proper "timing" for migration and evolution steps?

# Answer
> 3 votes
I assume you are simulating a null distribution. Are you investigation recombination??

My main advice is to use population genetic terminology rather than geomometry to describe your simulation (e.g. vertex is a migration event between allopatric populations). Migration is investigated between populations investigated via $F\_{ST}$. You appear to be assuming clonality, and population movement appear to follow a 'stepping stone model' (I might have misunderstand the question).

Mutation rate ... I'll have a guess at $10^5$ per site per year. Migration rate depends on the bacterial species being modeled, eg. MRSA spread in a very similiar fashion to that described in your model (along a moterway), in contrast soil bacteria don't adhere to allopatry at all (when I looked at them).

Fixing one parameter (population size) whilst investigation a second (migration) would be normal to assess the parameter space initially. Bacterial populations do undergo selective sweeps and using an appropriate growth model I would look at a standard population growth model (e.g. the one Bob May investigated in chaotic interactions within a differential equation). There are several population growth models and fixing migration whilst assessing the model would seem sensible. The parameter space in the combined model of migration, space and population size becomes complex.

---

Okay I get your model now. Its fine, the thing I'd be cautious about is the population density of the bacterial host because this would skew migration across the nodes of your graph. A quick literature review will demonstrate the importance, double paeroto log normal distribution and its cousin (the name escapes me .... ) are sort of key words that might flag up. Bacteria are not my active thing so I'm rusty.

---
Tags: phylogenetics, phylogeny, networks, bacteria, allele-frequency
---


# Question
Title: Biohackers Netflix - DNA to binary and video
I'm not sponsored or anything, just interested in their challenge to decipher their DNA code.

They encoded their first episode of "Biohackers" video/binary file to DNA code and said if we could decode it we can watch it (without Netflix). Here's their page: https://biohackersnetflix.com with description and download for the DNA sequence file. (Don't know if it's just in German or you can translate it. If questions regarding this page, ask me.)

The file is ~550MB small and contains 3.882.771 lines (not in fasta format). Every line has a length of 147 characters including primers at both ends (Illumina?). Here are the first 5 lines:

```
ACACGACGCTCTTCCGATCTCTCCCAGGGACAAAGGTTCTGCATTTGCAGCAAGACTCCTGTAGTGCTGCAGATTCTCTGGTTGGATAGTACGGCGTACATTTCTGTATTGTAGCACCATGGGGTAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTTAAGGCTTCGTAACAGATATTCTATATCGTCACATTGGTCTGAAGGAAGTCGCCTATAATCGCTCCTCTGTTTTTTAAAACTGCTATGGACCCGCTGTTCGGTGGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTCATGGTATAAGTGTTAAGGGTAATAACCACCTACCCCCCTCATTGCTCGTTTTTCCTGGAACCTTAACATTCGCAATAGCTAGCTGTTTCCTAGTAGAACCAAGGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTAGGATGTAGTCACAGGTCATTGTCATTAACTCAACCGAGGACATAACACTAAGTCCCACTAGGCCTGGATTCTCTAACGCGGTCTCTCTATTGGGGGAAGGGGTGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTTCTGGTAAGGCGGGTTGATATCAGTCACCTCCCTTTGAGCTAAAATACGATGGCGATTTAGTGTGAAACTAATAATGCTTGTCATACCAGCAGTACCGGATCGGGAGATCGGAAGAGCACACGTCT

```

I trimmed all the primers and tried to decode {A, C, G, T} considering every permutation {00, 01, 10, 11} as the obvious(?) decryption method (4! = 24 possible decodings) using python.

Then I hoped to get 1 of these 24 files loaded into VLC media player or something to be played, but it didn't work and every file seemed to be broken in the same way. I think I'm missing something here.

Can I assume that a text file containing only 0's and 1's should be playable in VLC if the DNA code is correctly decrypted?

(If I am wrong here, please tell or move me.)

//Edit: I converted all 24 files to ASCII to see if there's some kind of "video-like header". (All videos have some sort of description in their first lines if opened in text editor?) But there's just gibberish.

//Edit: I saw that every 84th sequence position has a "T", which is kind of weird. So I tried to run my script again with these T's removed, but still no solution.

//Edit: I searched for "AVI", "264", "codec" and some other strings in every video file I created and hexdumped. Nothing found. For clarification: I translated the DNA into every 24 binary and then into their ASCII representation following the 19 upvotes answer: https://stackoverflow.com/questions/7290943/write-a-string-of-1s-and-0s-to-a-binary-file. The 104 bases / 208 Bits (removed repetetive "T" and primers) are actually a multiple of 8 (respectively 26 Bytes) so I could be on the right way (even if not 32 Bytes?). De novo Assembly didn't work and I found no obvious ORF "genes" representing some kind of URL to the video or something which was a neat idea considering the video file would be only ~150MB. (See comments.)

# Answer
> 6 votes
I was quite curious when I saw that as well. I've spent more time than I'd care to admit trying the sort of things you have. I've gotten it decoded now, but I can't really claim any sort of victory, as the problem was pre-solved.

After struggling through some of the same experiments you did, I decided to take a closer look at their explainer video here: https://youtu.be/DMYgjOHgHxc

First I tried to decipher the blackboard sketch around 1:30, but it was quite vague. A closer investigation into the speaker, however, yielded some luck. Googling Dr. Reinhard Heckel brought me to his website here: http://www.reinhardheckel.com/ which shows his most recent publication – a paper on encoding digital data in DNA.

The encoding is relatively complex (something nearly impossible to stumble upon accidentally, I'd reckon) but fascinating. The fragments are indexed and have two, layered error-correcting codes. Perhaps making our job harder (but serving the practical purpose of minimizing homopolymers and unwanted annealing) is the fact that the data is XORed with a pseudo-random noise that shuffles the data. The paper obviously has all of the details.

From my skim, however, I stumbled across this Github repo linked in the paper: https://github.com/reinhardh/dna\_rs\_coding

As it happens, the README was just updated in this last week, describing how to decode the episode from the file provided. If you have Docker, you just need to copy-paste some commands.

The product is indeed the full and final episode – coming in at a dainty 63.1MB. How did they get 40 minutes of video that compressed? Well, in short, 720x360 at 24 fps. It looks a little awful and is in German only, but is certainly a cool little easter egg. If nothing else, I have a cool paper to read.

---
Tags: python, phylogenetics, encode
---


# Question
Title: 10X Illumina demultiplexing sample sheet issue
Also posted on biostars.

I am trying to use `cellranger` or `bcl2fastq` to convert the `.bcl` files that I got from single cell analysis run into `fastq` files for further analysis. I needed to generate `sample_sheet.csv` and so I used the following tool:

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/bcl2fastq-direct

which I copy-pasted into a text file using `vim`. `sample_sheet.csv` looks like that:

```
[Header]
EMFileVersion,4

[Reads]
26
8
98

[Data]
Lane,Sample_ID,Sample_Name,index,Sample_Project
1,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
1,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
1,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
1,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
1,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
1,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
1,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
1,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
2,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
2,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
2,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
2,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
2,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
2,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
2,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
2,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
3,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
3,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
3,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
3,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
3,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
3,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
3,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
3,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
4,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
4,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
4,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
4,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
4,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
4,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
4,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406

```

Then I ran `cellranger` using `bash` script:

```
#!/bin/bash

 FLOWCELL_DIR="/scratch/nv4e/kipnis/180403_NB501830_0158_AHN3LLBGX5"
 OUTPUT_DIR="/scratch/nv4e/kipnis/fastq"
 SAMPLE_SHEET_PATH="/scratch/nv4e/kipnis/sample_sheet_2.csv"

cellranger mkfastq --id="AHN3LLBGX5" --run=${FLOWCELL_DIR} --csv=${SAMPLE_SHEET_PATH} --output-dir=${OUTPUT_DIR}

```

It fails with the following output:

`_stderr` file shows:

and that my `sample_sheet.csv` can not be parsed:

It could be because `CLRF` needs to be converted into `LF`, but I tried two commands to do this unsuccessfully:

`tr -d '\r' < input > output` and `perl -pi -e 's/\r\n/\n/g' input`

From the following thread:

https://unix.stackexchange.com/questions/277217/how-to-install-dos2unix-on-linux-without-root-access

The above commands successfully strip the `\r` from the line endings, but the error from bcl2fastq remains.

What am I doing wrong? Any suggestions would be greatly appreciated.

> Update

I looked into the folder that is generated when the `cellranger` runs and found two `.csv` files generated by the `cellranger`. One of them, `samplesheet.csv`, contains these `^M` characters that are shown in the error:

> Update

Maybe one of the issues was that I was submitting regular `.sh` file on the `SLURM` cluster. However, I wrote the right `.slurm` file and submitted it yesterday with `sbatch` command:

However, today I got the message in my email:

And neither `_stderr`, nor `generate_fastq_id.err` or `generate_fastq_id.out` files contain anything. I thought initially that I set not enough time and it timed out, however, it is not true: I set the time to `24h` and it has been running for just `13h` before failing.

> Update

I searched over all of the files in the generated folder with the `find` linux command and actually found the error files that show the same error as before: issue with the `sample_sheet.csv` file. This error file is `AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/BCL2FASTQ_WITH_SAMPLESHEET/fork0/chnk0-u179acba3d8/_stderr`.

So, right now the issue seems to be the following: I am feeding in the right `sample_sheet.csv` without the `^M` symbols, but `cellranger` transforms it into one that contains `^M` and then uses it giving the error. Above I see also `which: no configureBclToFastq.pl` and `ERROR: bcl2fastq::common::Exception: 2018-Apr-09 13:33:14: Inappropriate ioctl for device`, so maybe my input parameters are just wrong to `cellranger` and it has nothing to do with the `sample_sheet.csv`?

> Update

I tried to use `--sample-sheet` instead of `--csv` since I am feeding in complex sample sheet. No luck.

> Update

Here you can find all the error document in text:

> which: no configureBclToFastq.pl in (/sfs/lustre/scratch/nv4e/cellranger/samtools\_new/1.6:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-fastq/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-ref/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/miniconda-cr-cs/4.3.21-miniconda-cr-cs-c9/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/bin:/sfs/lustre/scratch/nv4e/cellranger/lz4/v1.8.0:/sfs/lustre/scratch/nv4e/cellranger/martian-cs/2.3.2/bin:/sfs/lustre/scratch/nv4e/cellranger/STAR/5dda596:/sfs/lustre/scratch/nv4e/cellranger/samtools\_new/1.6:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-fastq/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-ref/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/miniconda-cr-cs/4.3.21-miniconda-cr-cs-c9/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/bin:/sfs/lustre/scratch/nv4e/cellranger/lz4/v1.8.0:/sfs/lustre/scratch/nv4e/cellranger/martian-cs/2.3.2/bin:/sfs/lustre/scratch/nv4e/cellranger/STAR/5dda596:/scratch/nv4e/cellranger:/scratch/nv4e/bcl2fastq/build/bin:/scratch/nv4e/spark/bin:/scratch/nv4e/scala/bin:/sfs/nfs/blue/nv4e/private/bin:/sfs/nfs/blue/nv4e/anaconda2/bin:/sfs/nfs/blue/nv4e/.local/bin:/usr/lib64/qt-3.3/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/opt/slurm/current/bin:/opt/slurm/current/sbin:/opt/singularity/current/bin:/opt/rci/bin:/opt/rci/sbin:/opt/nhc/current/sbin:/share/rci\_apps/common/bin:/share/resources/HPCtools/) BCL to FASTQ file converter bcl2fastq v2.20.0.422 Copyright (c) 2007-2017 Illumina, Inc. 2018-04-09 18:31:29 \[7fcdb39f37c0\] Command-line invocation: bcl2fastq --minimum-trimmed-read-length 8 --mask-short-adapter-reads 8 --create-fastq-for-index-reads --ignore-missing-positions --ignore-missing-filter --ignore-missing-bcls --use-bases-mask=Y26,I8,Y98 -R /scratch/nv4e/kipnis/180403\_NB501830\_0158\_AHN3LLBGX5 --output-dir=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE\_FASTQS\_CS/MAKE\_FASTQS/BCL2FASTQ\_WITH\_SAMPLESHEET/fork0/chnk0-u21f1cbe9bf/files/fastq\_path --interop-dir=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE\_FASTQS\_CS/MAKE\_FASTQS/BCL2FASTQ\_WITH\_SAMPLESHEET/fork0/chnk0-u21f1cbe9bf/files/interop\_path --sample-sheet=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE\_FASTQS\_CS/MAKE\_FASTQS/PREPARE\_SAMPLESHEET/fork0/chnk0-u09e7cbe891/files/samplesheet.csv -p 6 -r 6 -w 6 2018-04-09 18:31:29 \[7fcdb39f37c0\] INFO: Minimum log level: INFO 2018-04-09 18:31:29 \[7fcdb39f37c0\] INFO: Sample sheet: '/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE\_FASTQS\_CS/MAKE\_FASTQS/PREPARE\_SAMPLESHEET/fork0/chnk0-u09e7cbe891/files/samplesheet.csv' 2018-04-09 18:31:29 \[7fcdb39f37c0\] ERROR: bcl2fastq::common::Exception: 2018-Apr-09 18:31:29: Inappropriate ioctl for device (25): /scratch/nv4e/bcl2fastq/src/cxx/include/common/CsvGrammar.hpp(92): Throw in function bcl2fastq::common::CsvGrammarAttribute bcl2fastq::common::parseCsvData(Iterator, Iterator) \[with Iterator = \__gnu\_cxx::\__normal\_iterator \>; bcl2fastq::common::CsvGrammarAttribute = std::vector \> \>\] Dynamic exception type: boost::exception\_detail::clone\_impl std::exception::what: Could not parse the CSV stream text: ^M ^M \[Reads\]^M 26^M 98^M ^M ^M \[Data\]^M Lane,Sample\_ID,Sample\_Name,index,Sample\_Project,Original\_Sample\_ID^M 1,SI-GA-B4\_1,17R,ACTTCATA,Chromium\_20180409,SI-GA-B4\_1^M 1,SI-GA-B4\_2,17R,GAGATGAC,Chromium\_20180409,SI-GA-B4\_2^M 1,SI-GA-B4\_3,17R,TGCCGTGG,Chromium\_20180409,SI-GA-B4\_3^M 1,SI-GA-B4\_4,17R,CTAGACCT,Chromium\_20180409,SI-GA-B4\_4^M 1,SI-GA-B5\_1,19RL,AATAATGG,Chromium\_20180409,SI-GA-B5\_1^M 1,SI-GA-B5\_2,19RL,CCAGGGCA,Chromium\_20180409,SI-GA-B5\_2^M 1,SI-GA-B5\_3,19RL,TGCCTCAT,Chromium\_20180409,SI-GA-B5\_3^M 1,SI-GA-B5\_4,19RL,GTGTCATC,Chromium\_20180409,SI-GA-B5\_4^M 2,SI-GA-B4\_1,17R,ACTTCATA,Chromium\_20180409,SI-GA-B4\_1^M 2,SI-GA-B4\_2,17R,GAGATGAC,Chromium\_20180409,SI-GA-B4\_2^M 2,SI-GA-B4\_3,17R,TGCCGTGG,Chromium\_20180409,SI-GA-B4\_3^M 2,SI-GA-B4\_4,17R,CTAGACCT,Chromium\_20180409,SI-GA-B4\_4^M 2,SI-GA-B5\_1,19RL,AATAATGG,Chromium\_20180409,SI-GA-B5\_1^M 2,SI-GA-B5\_2,19RL,CCAGGGCA,Chromium\_20180409,SI-GA-B5\_2^M 2,SI-GA-B5\_3,19RL,TGCCTCAT,Chromium\_20180409,SI-GA-B5\_3^M 2,SI-GA-B5\_4,19RL,GTGTCATC,Chromium\_20180409,SI-GA-B5\_4^M 3,SI-GA-B4\_1,17R,ACTTCATA,Chromium\_20180409,SI-GA-B4\_1^M 3,SI-GA-B4\_2,17R,GAGATGAC,Chromium\_20180409,SI-GA-B4\_2^M 3,SI-GA-B4\_3,17R,TGCCGTGG,Chromium\_20180409,SI-GA-B4\_3^M 3,SI-GA-B4\_4,17R,CTAGACCT,Chromium\_20180409,SI-GA-B4\_4^M 3,SI-GA-B5\_1,19RL,AATAATGG,Chromium\_20180409,SI-GA-B5\_1^M 3,SI-GA-B5\_2,19RL,CCAGGGCA,Chromium\_20180409,SI-GA-B5\_2^M 3,SI-GA-B5\_3,19RL,TGCCTCAT,Chromium\_20180409,SI-GA-B5\_3^M 3,SI-GA-B5\_4,19RL,GTGTCATC,Chromium\_20180409,SI-GA-B5\_4^M 4,SI-GA-B4\_1,17R,ACTTCATA,Chromium\_20180409,SI-GA-B4\_1^M 4,SI-GA-B4\_2,17R,GAGATGAC,Chromium\_20180409,SI-GA-B4\_2^M 4,SI-GA-B4\_3,17R,TGCCGTGG,Chromium\_20180409,SI-GA-B4\_3^M 4,SI-GA-B4\_4,17R,CTAGACCT,Chromium\_20180409,SI-GA-B4\_4^M 4,SI-GA-B5\_1,19RL,AATAATGG,Chromium\_20180409,SI-GA-B5\_1^M 4,SI-GA-B5\_2,19RL,CCAGGGCA,Chromium\_20180409,SI-GA-B5\_2^M 4,SI-GA-B5\_3,19RL,TGCCTCAT,Chromium\_20180409,SI-GA-B5\_3^M 4,SI-GA-B5\_4,19RL,GTGTCATC,Chromium\_20180409,SI-GA-B5\_4^M

# Answer
So, we figured it out. If you delete everything until the `[Data]` tag in the sample sheet, leaving the `[Data]` tag, it starts working. Here is the full `sample_sheet.csv` that is working:

```
[Data]
Lane,Sample_ID,Sample_Name,index,Sample_Project
1,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
1,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
1,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
1,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
1,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
1,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
1,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
1,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
2,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
2,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
2,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
2,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
2,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
2,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
2,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
2,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
3,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
3,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
3,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
3,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
3,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
3,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
3,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
3,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
4,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
4,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
4,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
4,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
4,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
4,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
4,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
4,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409

```

It is very unfortunate and disappointing that `10X` website is generating wrong sample sheet. Still we got very low mapping rates and trying to figure out what went wrong.

> Update

I was supplying `--sample-sheet=` property, but it should be `--samplesheet=`. Now it is running fine (except for very low mapping rates) with all the headers on.

> Update

It turned out that we were given the wrong indexes, so we made a sample sheet using sample sheet generator again and even with `--samplesheet=` correct parameter it was not running. Convertion of `CLRF` with `perl` or `tr` did not help. However, concise format of the sample sheet worked out well which looks like the following:

```
Lane,Sample,Index
1-4,17R,SI-GA-B8
1-4,19RL,SI-GA-B9

```

It needs to be supplied with `--csv=` parameter.
> 2 votes

# Answer
Apart from supporting only Windows-style and Linux-style line endings and a restricted character set as mentioned above, the sample sheet has to be in the correct encoding which is UTF-8 (or its subset US-ASCII). Many command lines support these lookup and conversion commands:

```
# lookup
file -i my_sample_sheet.csv

# conversion
iconv -f ISO-8859-15 -t UTF-8 in.csv > out.csv

```
> 0 votes

---
Tags: scrnaseq, 10x-genomics
---


title: Lifting from hg38 to hg37

text:
# Question
Title: Lifting from hg38 to hg37
I am having strange issues lifting this dataset to hg37 from hg38. I am using the HGDP WGS data from ftp://ngs.sanger.ac.uk/production/hgdp/hgdp\_wgs.20190516/. The data set split by chromosome, so I first converted them to BED then merged them to a single file using plink1.9. However, when trying to lift to hg37/19 using either liftOver or CrossMap v0.5.1 I get strange errors. I tried both the chain from ensemble.org and from UCSC, I tried them both compressed and uncompressed. I am starting to think I need to try the conversion on the VCF files but these are so large I would prefer not to duplicate them unnecessarily.

# Data

The data looks fine. If we look at the BIM:

```
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38_chrpos.bim
1   1:10153 0   10153   G   A
1   1:10163 0   10163   C   T
1   1:10180 0   10180   C   T
1   1:10250 0   10250   C   A
1   1:10257 0   10257   C   A
1   1:10291 0   10291   T   C
1   1:10297 0   10297   T   C
1   1:10327 0   10327   C   T
1   1:10330 0   10330   A   C
1   1:10333 0   10333   T   C
user@desktop:/media/luks8tb1/data/genomics/HGDP$ tail hgdp_hg38_chrpos.bim
Y   Y:56887228  0   56887228    C   T
Y   Y:56887285  0   56887285    C   A
Y   Y:56887316  0   56887316    A   G
Y   Y:56887463  0   56887463    A   C
Y   Y:56887491  0   56887491    A   C
Y   Y:56887503  0   56887503    T   G
Y   Y:56887583  0   56887583    G   A
Y   Y:56887631  0   56887631    T   C
Y   Y:56887837  0   56887837    A   G
Y   Y:56887844  0   56887844    C   T

```

The data are in chrpos variant names here but I had the same result using the original rsid.

# Chains

Chains look fine, and I am able to read them into R without problems:

```
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head GRCh38_to_GRCh37.chain
chain 1 1 248956422 + 10000 297968 1 249250621 + 10000 267719 2
167417 80249 50000
40302

chain 1 1 248956422 + 347968 501617 1 249250621 - 248779253 248932902 3
153649

chain 1 1 248956422 + 585988 12949384 1 249250621 + 521368 13009210 4
1044707 1 0
3709 3 0
user@desktop:/media/luks8tb1/data/genomics/HGDP$ tail GRCh38_to_GRCh37.chain

chain 1 chrY 57227415 + 56821509 57217415 Y 59373566 + 58967656 59363566 11183
85168 0 1
78580 0 1
63496 0 1
3158 1 0
30382 0 1
54147 0 1
80974

```

# CrossMap

```
user@desktop:/media/luks8tb1/data/genomics/HGDP$ CrossMap.py bed GRCh38_to_GRCh37.chain hgdp_hg38_chrpos.bed hgdp_hg37_chrpos.bed
@ 2020-08-25 20:04:26: Read the chain file:  GRCh38_to_GRCh37.chain
Traceback (most recent call last):
  File "/home/user/.local/bin/CrossMap.py", line 166, in <module>
    crossmap_bed_file(mapTree, in_file, out_file)
  File "/home/user/.local/lib/python3.6/site-packages/cmmodule/mapbed.py", line 32, in crossmap_bed_file
    for line in ireader.reader(inbed):
  File "/home/user/.local/lib/python3.6/site-packages/cmmodule/ireader.py", line 27, in reader
    yield l.decode('utf8').strip().replace("\r", "")
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xdd in position 20: invalid continuation byte

```

# liftOver

```
user@desktop:/media/luks8tb1/data/genomics/HGDP$ liftOver hgdp_hg38.bed hg38ToHg19.over.chain hgdp_hg37.bed hgdp_hg37_errors.bed
Reading liftover chains
Mapping coordinates
Data format error: expecting at least 3 fields in BED file (hgdp_hg38.bed)

```

# Trying the PED file with liftOver

I converted BED file to PED with plink1.9:

```
plink --bfile hgdp_hg38 --recode --out hgdp_hg38

```

Then:

```
user@desktop:/media/luks8tb1/data/genomics/HGDP$ liftOver hgdp_hg38.ped hg38ToHg19.over.chain hgdp_hg37.ped hgdp_hg37_errors.ped
Reading liftover chains
Mapping coordinates
invalid unsigned integer: "HGDP00001"

```

There's a few question on Google about this, but nothing that looks useful.

CrossMap does not support PED.

# Answer
BED format can refer to the ".txt" file format that stores chromosomes and positions of genomic features as well as some other things like name. `bed` format can also refer to the plink format, which belongs to bed, bim, fam trio. The plink bed format stores all the genetic variants, but the program needs to refer to bim file for their positions.

Genome liftover obviously only needs the coordination of each SNPs, and it does not need individual genotypes at all. In fact, the bed format in PLINK does not provide any useful information. Knowing this means that providing the plink bed file is wrong. And this is pretty clear from the error message.

Moreover, since the three plink file always needs to be used together, most programs refer to them as **plink file format**. So whenever people say BED format, it is the txt file format.
> 2 votes

# Answer
Not sure about crossMap but liftOver requires the BED format where column 1 is the chromosome, column 2 the 0-based start and column 3 the end coordinate. Additional columns are optional and can be used to store the other columns you have. Simply convert BIM to BED and run liftOver.
> 0 votes

---
Tags: plink, liftover
---


# Question
Title: Error When Using biocLite as an installer in rpy2 python library
I am running into an error when I tried to run the following lines of code, the `error` at line 3:

```
1 from rpy2.robjects.packages import importr
2 base = importr('base')
3 base.source("http://www.bioconductor.org/biocLite.R")
4 biocinstaller = importr("BiocInstaller")
5 biocinstaller.biocLite("seqLogo")

```

`rpy2.rinterface_lib.embedded.RRuntimeError: Error: With R version 3.5 or greater, install Bioconductor packages using BiocManager; see https://bioconductor.org/install`

How can I use `BiocManager` installer instead of `biocLite`  and reflect that on my code? I couldn't find any documentation or source that talk about this error from within Python.

Thank you in advance

# Answer
> 2 votes
You're using an outdated version of R. You can use the R 4.0 from `conda-forge` channel, just create an environment with `r-base` and then install `rpy2`, the latest version compatible with R 4.0 is available through `pip`

```
conda create -n rpy2 -c conda-forge r-base
conda activate rpy2
pip install rpy2

```

Now you can use rpy2 with R 4.0 in your python code.

```
import rpy2.robjects as robjects
robjects.r('install.packages("BiocManager", '
           'repos="http://cran.r-project.org")')
robjects.r('BiocManager::install("seqLogo")')
robjects.r('library(seqLogo)')

```

Remember deactivate your env with `conda deactivate`.

---
Tags: r, python, rpy2
---


# Question
Title: Microarray Meta-Analysis and Cross-Platform
I have a microarray dataset with Agilent platform, all of the samples are cancerous, I need normal or control samples to compare with cancerous samples, can I combine control samples from another platform (for example Illumina) to Agilent platform? and then remove the batch effect (by `ComBat` function in the sva package).

I would be grateful if you could guide me to do it.

# Answer
> 4 votes
In short the answer is no because you have to standardise between the arrays. Thus both arrays would need to run exactly the same control sample and you would then perform statistical standardisation (its a formal term), which is a varience measure where you look at the variation around the 'mean' for all samples in the array and here the 'mean' is the standard sample. If that can be done for both arrays at that point they can be combined.

The problem is intensity per run varies and you have to control against this

---

In this scenario you can certainly combine samples by assessing whether the control samples between runs approximates to identity. What you could use is a uniform distribution to assess this, variance must MUCH smaller than the mean for the STANDARD sample(s). I don't really do uniform distribution statistics, so I'd need to research this, because 'uniform' means identity to standard deviations against the uniform distribution is a bit meaningless. However again the uniform distribution is denoted by a very variance in comparison to the mean. The statistic of dividing mean by standard deviation/variance from memory is called the Index of Association (I could be wrong). If this value is at an extreme you're good to go \<\<\<1, if it is remotely close to `1` you're in trouble, if it is \>1 you need to use standardisation.

Presumably you have several standards to ensure a mean is meaningful. A T-test and non-parametric equivalents is not stringent enough. If there is a difference assess a correction factor and rerun the uniform distribution check. The same correction will be applied to all samples.

The undisputed method to combine is standardisation, but accept in this scenario the differences between runs is potentially very small, so here what you are attempting is essentially to combine raw data or at least keep one of the runs as pure raw data.

---
Tags: r
---


# Question
Title: Sort reads in BAM file based on presence of specific deletion?
I have an indexed BAM file containing long-read sequencing data and I'd like to split the reads contained within into those with a known deletion and those without the deletion (I have the deletion coordinates available to me) when mapped against the hg38 genome. Ideally I'd like to end up with two BAM files (one with reads containing the deletion and one without) for downstream analysis. I've searched a bit online and haven't really found any tools that do exactly that, but this seems like a common enough task that there has to be something.

Thanks in advance for any help/guidance!

**UPDATE**

In addition to @gringer's great answer, I decided to write a small python script myself that solves my issue using pysam, the fruits of which can be found here. If anyone has any feedback or suggestions on it I'll happily take them, as I'm still very new to this. It's undergone very little testing, so user be warned.

# Answer
> 1 votes
I was asked to write some code sometime near the end of last year to do this, so that viral sequences could be split based on variants at a particular location. The code I wrote doesn't actually do the splitting - it creates SAM read groups based on those variants - but that splitting can be done after processing with this script using `samtools split`.

As with most of my scripts, it hasn't had much testing / use:

https://gitlab.com/gringer/bioinfscripts/-/blob/master/samVarSplitter.pl

usage:

```
samtools view -h mapped_reads.bam | ./samVarSplitter.pl [-ref <ref>] [-pos <int>] [options]

```

example:

```
samtools view -h mapped_reads.bam | ./samVarSplitter.pl -ref Wuhan-Hu-1 -pos 11083

```

If the input file contains a lot of reads that don't include the variant location, I'd recommend pre-filtering reads to only include those that cover that location, because expanding the CIGAR string and storing a record of the positions can take some time.

Note: currently the position argument only supports a single location. I could extend this to multiple positions (i.e. a position range) if that were a desired feature.

# Answer
> 2 votes
The CIGAR string should tell you if a read has a deletion in it. So filter the bam files based on that.

---
Tags: sam, long-reads, indel
---


title: Dealing with duplicate keys

text:
# Question
Title: Dealing with duplicate keys
I got the following error:

> attempting to use to\_dict to index a multiple sequence FASTA file with duplicate keys.

Are there methods for working around duplicate keys in biopython? I was thinking to simply remove the duplicate entries, but I would need to delete the entire entry associated with the keys.

# Answer
Yes, Biopython deliberately does not allow duplicate keys with `.to_dict()`, e.g.

Input file:

```
$ cat test.fa
>1
ATG
>2
AAA
>1
CG

```

Code:

```
from Bio import SeqIO

print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta")))

```

Traceback:

```
  File "script.py", line 3, in <module>
    print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta")))
  File "/path/lib/python3.7/site-packages/Bio/SeqIO/__init__.py", line 789, in to_dict
    raise ValueError("Duplicate key '%s'" % key)
ValueError: Duplicate key '1'

```

However, you can easily write your own function, e.g.:

```
from Bio import SeqIO

def to_dict_remove_dups(sequences):
    return {record.id: record for record in sequences}

print(to_dict_remove_dups(SeqIO.parse("test.fa", "fasta")))

```

Output:

```
{'1': SeqRecord(seq=Seq('CG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[]), '2': SeqRecord(seq=Seq('AAA', SingleLetterAlphabet()), id='2', name='2', description='2', dbxrefs=[])}

```

Alternatively, you Biopython's `.to_dict()` accepts a `key_function` argument, e.g. if you wish to make all the headers unique

```
from Bio import SeqIO
from itertools import count

c = count()

print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta"),
                    key_function = lambda x: f"{x.id}_seqNum{next(c)}"))

```

Output:

```
{'1_seqNum0': SeqRecord(seq=Seq('ATG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[]), '2_seqNum1': SeqRecord(seq=Seq('AAA', SingleLetterAlphabet()), id='2', name='2', description='2', dbxrefs=[]), '1_seqNum2': SeqRecord(seq=Seq('CG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[])}

```
> 2 votes

---
Tags: python
---


# Question
Title: COVID-19 GWAS: what is a cross-replicating association?
In a genomic study of patients infected with SARS-CoV-2, the authors detected cross-replicating associations with rs11385942 at locus 3p21.31 , the association signal spanned the genes SLC6A20, LZTFL1, CCR9, FYCO1, CXCR6 and XCR1

*I would like to know what's a cross-replication association.*

Source pdf is here

# Answer
Replication in GWAS is an analytical term, approximating to 'robustness'. It is not biological phenomenon of a genome replicating. What the authors are saying is these genes associate with an allelic variant(s) which correlate with the disease phenotypic of COVID-19.

Three of those genes are very interesting indeed notably CCR9, CXCR6 and XCR1, because these are chemokine receptors. The later for example is associated with regulating between the innate and adaptive (CD8+, i.e. cytotoxic T-cell) immune responses (described in the image below). The innate/adaptive interface is an extremely trendy area in immunology.

The authors are implying the immunological signalling associated with one or more of these chemokines is fundamental to combating COVID-19 infection. 'Replication' is simply the GWA geek for 'its real'.
> 2 votes

---
Tags: genomics, gene-expression, covid-19, sars-cov-2, coronavirus
---


title: What are phantom peaks in ChIP-seq?

text:
# Question
Title: What are phantom peaks in ChIP-seq?
I've seen two uses of this term that seem to refer to completely different phenomena. Is my understanding correct?

1. In Jain et al. "Active promoters give rise to false positive ‘Phantom Peaks’ in ChIP-seq experiments" (PMC4538825), a phantom peak is a region on the genome where lots of reads align in a control experiment where there is none of the protein being IP'd, but not in a control experiment without any antibody. The authors (speculatively) attribute this to general stickiness of transcription factors and the transcriptional machinery, which can bind with potentially any antibody and get IP'd.
2. The ENCODE quality metrics (link) say "a cross-correlation phantom-peak is ... observed at a strand-shift equal to the read length." They are referring to a high correlation of x\[i\] and y\[i + r\] where r is the read length, x\[i\] contains counts of tags aligning to the Watson strand at coordinate i, and y\[i+r\] contains counts of tags aligning to the Crick strand at coordinate i+r. They attribute this phenomenon to variations in mappability. (I don't understand their explanation at all, but that's a separate question.)

# Answer
> 1 votes
This is almost right, but it's missing a little context. In ChIP-seq, it's typical to shift the positive and negative reads because otherwise the peaks end up offset from each other. The amount of shift is determined by trying different shifts and measuring how well the read densities correlate between positive and negative strands. In the question, the first sense talks about a peak on the plot of read density vs genomic position. The second has different x and y axes because it's talking about the cross-correlation between strands. X is the offset between positive and negative strands, and y is the correlation observed at that offset. Look at figures 4D,E in this ENCODE paper for visuals.

Landt, S. G., Marinov, G. K., Kundaje, A., Kheradpour, P., Pauli, F., Batzoglou, S., ... & Chen, Y. (2012). ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome research, 22(9), 1813-1831. https://genome.cshlp.org/content/22/9/1813.full.pdf

---
Tags: chip-seq
---


# Question
Title: How to use pyopenms.FeatureFindingMetabo() for peak picking in pyOpenMS?
I found this documentation on pyOpenMS mass spec library. I want to try different peak picking algorithms. Unfortunately, the docs are not very detailed and don't show how to use the other peak picking algorithms and the codebase itself is also quite cryptic for me.

Any idea how to use, for example, the FeatureFindingMetabo algorithm insead of the centroided one?

```
from urllib.request import urlretrieve
# from urllib import urlretrieve  # use this code for Python 2.x
gh = "https://raw.githubusercontent.com/OpenMS/OpenMS/develop"
urlretrieve (gh +"/src/tests/topp/FeatureFinderCentroided_1_input.mzML", "feature_test.mzML")

from pyopenms import *

# Prepare data loading (save memory by only
# loading MS1 spectra into memory)
options = PeakFileOptions()
options.setMSLevels([1])
fh = MzMLFile()
fh.setOptions(options)

# Load data
input_map = MSExperiment()
fh.load("feature_test.mzML", input_map)
input_map.updateRanges()

ff = FeatureFinder()
ff.setLogType(LogType.CMD)

# Run the feature finder
name = "centroided"
features = FeatureMap()
seeds = FeatureMap()
params = FeatureFinder().getParameters(name)
ff.run(name, input_map, features, params, seeds)

features.setUniqueIds()
fh = FeatureXMLFile()
fh.store("output.featureXML", features)
print("Found", features.size(), "features")

```

# Answer
> 1 votes
This is a draft answer, if the OP still need needs I will extend it later.

Me too I found the documentation and examples of `pyopenms` very scarce. At the same time OpenMS is huge, sometimes you find key information not in various other sources instead of the official documentation. This is especially true if it comes to metabolomics.

We built a Python OpenMS lipidomics LC MS/MS preprocessing pipeline, contained in this file: https://github.com/saezlab/lipyd/blob/master/src/lipyd/msproc.py

The top class is the one called `MSPreprocess`, it implements the workflow, the other classes implement its steps.

Briefly, you start from profile mode mzML files, then

1. Peak picking with PeakPickerHiRes results centroided data
2. FeatureFindingMetabo assembles the traces in a metabolomics compatible way
3. A MapAlignment algorithm alignes the traces across scans
4. A FeatureGrouping algorithm alignes the features across multiple experiments

Our software mentioned above (`lipyd`) outputs data frames from the m/z, RT and intensity.

About `OpenMS` I can recommend to read the `C++` code and docs. The Python classes and methods are automatically built minimal wrappers, everything works the same way as in `C++.`

# Answer
> 0 votes
My minimal procedure based on the code that @deeenes shared for a single file. Comments to improve this solution are highly appreciated:

```
files = list_of_filenames

feature_map = oms.FeatureMap()
mass_traces = []
mass_traces_split = []
mass_traces_filtered  = []

exp = oms.MSExperiment()

options = oms.PeakFileOptions()
options.setMSLevels([1])

fh = oms.MzXMLFile()
fh.setOptions(options)

# Peak map
peak_map = oms.PeakMap()

file = files[0]

print('# Filename:', basename(file))
fh.load(file, exp)

print('# Spectra:', len( exp.getSpectra() ))
print('# Chromatograms:', len( exp.getChromatograms() ) )

for chrom in exp.getChromatograms():
    peak_map.addChrom(chrom)

for spec in exp.getSpectra():
    peak_map.addSpectrum(spec)

mass_trace_detect = oms.MassTraceDetection()
mass_trace_detect.run(peak_map, mass_traces, 100000)

print('# Mass traces:', len(mass_traces) )

elution_peak_detection = oms.ElutionPeakDetection()
elution_peak_detection.detectPeaks(mass_traces, mass_traces_split)
print('# Mass traces split:', len(mass_traces_split) )

feature_finding_metabo = oms.FeatureFindingMetabo()
feature_finding_metabo.run(
            mass_traces_split,
            feature_map,
            mass_traces_filtered)

print('# Mass traces filtered:', len(mass_traces_filtered) )
print('# Features:', feature_map.size() )

oms.FeatureXMLFile().store('FeatureFindingMetabo.featureXML', feature_map)

>>>
# Filename: CA_1.mzXML
# Spectra: 989
# Chromatograms: 0
# Mass traces: 18476
# Mass traces split: 25194
# Mass traces filtered: 0
# Features: 23549

```

---
Tags: python
---


# Question
Title: How are different forcefields developed/evaluated?
The Rosetta software can use a variety of different energy functions:

* Ref2015
* GenPot
* Beta-Nov16

How are these energy functions developed and validated? Is it common to do a single protein design with several of these energy functions?

# Answer
> 1 votes
Some force-fields rely solely on theoretical physical principles, while others are tweaked specifically for the goal of accurately depicting certain structures like sheets and helices or movement/diffusion aspects at different theoretical temperatures of a set of well-studied proteins. (Top-Down vs. Bottom-Up)

I haven't done MD myself but I am in an MD group and people seem to normally stick to a single force-field for a given MD project but it happens that they do very short simulation runs with force-fields just as a test. They definitely don't plan on doing the full simulation on different force-fields unless the comparison of those is their main goal.

# Answer
> 1 votes
The Rosetta ref2015 is well described in two papers, which I wish I had read much earlier, see this JChIM paper.

It's a statistical forcefield, where the terms are weighted to best match the empirical total energy in kcal/mol and the paper does a swell job at explaining what they do.

For several Rosetta experiments, you switch between *full-atom* and *centroid* (=coarse grain). These have different weights, as you can see from the files in the database folder or via pyrosetta, but majorly they have different "params" (topologies) for the residues as the centroid lacks many atoms as it's a simplification. In some instances, you may switch between *dihedral* (normal) to *cartesian* (explained in the FastRelax paper), wherein the former used internal coordinates (e.g. atom A is 1.5 Å away from B and A,B,C is a 107ª angle etc —cf. a params file `ICOOR` block), while the later is cartesian coordinates ( A is at \[x, y, z\] —cf. PDB file).

The GenPot uses its own `AtomTypes`, which means it uses its own params for each residue —like centroid/fa. The topology of a residue (ResidueType) is based on its conformation (ICOOR), bonding and atom definition, which is atom name, gasteiger-marsili charge and AtomType (an amalgamation of element, sp-hybridization, size, H-bond donor/acceptor etc.) These AtomTypes are better for ligands as many groups behave different from peptide atoms —say, a nitrile nitrogen is definitely different from a deprotonated aromatic nitrogen.

Some application, such as hydration, ligand docking etc. have specifically tweaked scorefunctions for use by their movers. FastRelax is a nice example as the cycle weighting changes are spat out by the tracer at vertiginous speed, where the protocol is:

```
repeat 5
coord_cst_weight 1.0
scale:fa_rep 0.040
repack
scale:fa_rep 0.051
min 0.01
coord_cst_weight 0.5
scale:fa_rep 0.265
repack
scale:fa_rep 0.280
min 0.01
coord_cst_weight 0.0
scale:fa_rep 0.559
repack
scale:fa_rep 0.581
min 0.01
coord_cst_weight 0.0
scale:fa_rep 1
repack
min 0.00001
accept_to_best
endrepeat

```

As you can see, the repulsion term of the Lenard-Jones is going from 0.04 up. These values were chosen not for thermodynamic reasons, but worked best for the application (energy minimisation of the structure).

However, this switching is all internal. Scores are always reported with the same `scorefunction` and comparing between different ones is not advised even for the minority of `scorefunction`s calibrated to be in kcal/mol.

---
Tags: protein-structure, computational-biochemistry, rosetta
---


# Question
Title: FPKM, FPKM-UQ, TPM or counts: How do I know which kind of unit should I use?
I'm trying to download data from the TCGA for gene expression analyses in R, but I'm in doubt if I should use FPKM, FPKM-UQ or counts? When the dataset is in counts, I suppose it's raw data, isn't it? So what's the best unit to compare multiple datasets?

I'm planning to use limma or Dseq2 for GE analysis and found that with Dseq2 I need to use count(non-normalized???) data... is that correct? so what's the best package and working strategy?

thank you very much, Fabiano

# Answer
> 3 votes
Type of data you need depends on the downstream applications and since you would like to carry out DEA with `DESeq2`, you would need raw counts (non-normalized).

There are many ways to import/download TCGA data, one such tool, the `TCGAbiolinks` package gives a nice interface for not only downloading the read count (or pre-processed) data but also associated clinical data.

# Answer
> 1 votes
You should use "counts" from the list of gene expression metrics that you provided rather than either FPKM or FPKM-UQ. FPKM is a normalized gene expression metric and an acronym for fragments per kilobase per million mapped reads. FPKM normalizes for both sequencing depth and genome size. Those specific packages want non-normalized counts as the expression metric. This means the number of reads mapped to a particular gene or feature. This is a slightly ambiguous metric but it can be the number of single-end reads mapped to a loci, the amount of pair end reads counted towards a loci, or the amount of super-reads (the theoretical sequence that spans the pair of paired-end reads which includes the presumed sequence between the first and second strand of a set of pair end reads) counted toward a loci. These values should be integers based on the assumptions of RNA-Seq and microarray technology which measures RNA transcripts present in a sample and precludes transcripts from spanning multiple genes. It does not matter which of these metrics you use as long as you are consistent. Those Bioconductor packages provide several different normalization methods but rely on read counts as their input.

---
Tags: r, gene-expression
---


# Question
Title: wtdbg2: practical implications of k-mer fsize and psize choice
I am using wtdbg2 2.3 to assemble a human genome (sequenced on PromethION from a cell line). I filtered out reads with low average quality, and now I am trying to determine the parameters that will optimize overall assembly quality.

The wtdbg2 assembler has two k-mer related settings, documented as follows.

```
 -k <int>    Kmer fsize, 0 <= k <= 25, [0]
 -p <int>    Kmer psize, 0 <= p <= 25, [21]
             k + p <= 25, seed is <k-mer>+<p-homopolymer-compressed>

```

The requirement that k + p is \<= 25\** suggests that these parameters can be set simultaneously. However, in all of the presets for the various long read sequencing technologies, one of the two settings is always set to 0.

My questions are as follows.

1. Are `k` and `p` ever intended to be used simultaneously?
2. If so, how does the behavior change with `-k 2 -p 19` versus `-k 19 -p 2`?
3. I'm familiar with some common implications of choosing k size for k-mer based analyses of Illumina reads. (Smaller k-sizes are more sensitive but also more subject to spurious matches. Larger k-sizes are more specific but also more likely to include a sequencing error.) Do similar implications apply for `-k` and `-p` with wtdbg2?
4. BONUS: Is there an explanation anywhere of how `-k` and `-p` are used? It looks to me like `-k` is related to k-mers in the traditional sense, and `-p` is a homopolymer compressed k-mer. Are seeds made by concatenating the regular k-mer and the compressed k-mer? It's possible I missed some important details in the preprint but I haven't been able to make the connection yet.

---

\**Incidentally, the documentation here is wrong. During runtime wtdbg2 enforces a constraint of \<= 23, not 25. This matches the description in the bioRxiv preprint stating that 46 bits (2 bits per base) are used to store each k-mer.

# Answer
> 1 votes
There is some description of homocompressed polymers (-p kmer) in this paper: https://arxiv.org/pdf/1708.01492.pdf

---
Tags: genome, nanopore, assembly, long-reads
---