text:
How to use the contrast function in ebayes for making different comparisons
The design matrix i have is this I would like to know how to use they way its done in deseq2 where we can use the contrast function to make particular comparison
The code im running to make differential testing is this
tcgaGlm = lmFit(as.matrix(tcgaExpr), design = tcgaDesign)
tcgaGlm = eBayes(tcgaGlm)

Now I will show the part of the data frame and the question Im trying to ask
head(tcgaDesign)
             Offset CEBPA DNMT3A FAM5C FLT3_ITD FLT3_TKD IDH1 IDH2 KIT KRAS NPM1 NRAS PHF6 PTPN11 RAD21 RUNX1 SMC1A SMC3 STAG2 TET2 TP53 U2AF1 WT1
TCGA-AB-2803      1     0      0     0        0        0    0    0   0    0    0    0    0      0     0     0     0    0     0    0    0     0   0
TCGA-AB-2805      1     0      0     0        0        0    0    1   0    0    0    0    0      0     0     1     0    0     0    0    0     0   0
TCGA-AB-2806      1     0      0     0        0        0    0    0   0    0    0    0    0      0     0     0     0    0     0    0    0     0   0
TCGA-AB-2807      1     0      0     0        0        0    0    1   0    0    0    0    0      0     0     1     0    0     0    0    0     0   0
TCGA-AB-2808      1     1      0     0        0        0    0    0   0    0    0    1    0      0     0     0     0    0     0    0    0     0   0
TCGA-AB-2810      1     0      0     0        0        0    0    1   0    0    1    0    0      0     0     0     0    0     0    0    0     0   0
               complex  minus5_5q   minus7q  plus8_8q     plus21    t_15_17     t_8_21 inv16_t16_16
TCGA-AB-2803 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000 1.00000000 0.00000000   0.00000000
TCGA-AB-2805 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000   0.00000000
TCGA-AB-2806 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000 0.00000000 1.00000000   0.00000000
TCGA-AB-2807 0.1410256 0.08974359 0.1217949 0.1153846 1.00000000 0.09615385 0.04487179   0.04487179
TCGA-AB-2808 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000   0.00000000
TCGA-AB-2810 0.1410256 0.08974359 0.1217949 0.1153846 0.05063291 0.09615385 0.04487179   0.04487179

Here if I have to make a comparison between CEBPA vs DNMT3A mutation in samples which carries CEBPA mutation and other which carries DNMT3A what or how should I do it.
To begin with I followed this manual I m still not sure once you run the model do i need to form a design matrix which comprises of sample which contain mutation in CEBPA and DNMT3A ? or I can go ahead with this experimental design
dput(tcgaDesign)
structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,
1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1,
1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0,
1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,
0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0,
0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0,
0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,
1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1,
0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1,
1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0.141025641025641,
0, 0.141025641025641, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0.141025641025641, 0, 0, 0, 0.141025641025641, 0, 0,
0, 1, 0, 0, 0, 0.141025641025641, 0.141025641025641, 0.141025641025641,
0, 0, 0, 0, 0.141025641025641, 0, 0, 0.141025641025641, 1, 0,
1, 0, 0, 1, 0, 0, 0, 0.141025641025641, 0, 0, 1, 0, 0, 0, 1,
0, 0, 0, 0.141025641025641, 0, 0, 0, 0, 0, 1, 0, 0, 0.141025641025641,
0, 0, 1, 0.141025641025641, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0,
1, 0, 0, 0, 1, 0.141025641025641, 1, 0, 0, 1, 0, 0, 0, 0, 0,
1, 0, 0.141025641025641, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1,
0.141025641025641, 0.141025641025641, 0, 1, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0.0897435897435897, 0, 0.0897435897435897, 0, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0897435897435897,
0, 0, 0, 0.0897435897435897, 0, 0, 0, 1, 0, 0, 0, 0.0897435897435897,
0.0897435897435897, 0.0897435897435897, 0, 0, 0, 0, 0.0897435897435897,
0, 0, 0.0897435897435897, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0.0897435897435897,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.0897435897435897, 0, 0, 0, 0,
0, 1, 0, 0, 0.0897435897435897, 0, 0, 0, 0.0897435897435897,
0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0.0897435897435897,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0897435897435897, 0, 0, 0,
1, 0, 0, 0, 1, 0, 1, 0, 1, 0.0897435897435897, 0.0897435897435897,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.121794871794872, 0, 0.121794871794872,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.121794871794872,
0, 0, 0, 0.121794871794872, 0, 0, 0, 1, 0, 0, 0, 0.121794871794872,
0.121794871794872, 0.121794871794872, 0, 0, 1, 0, 0.121794871794872,
0, 0, 0.121794871794872, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0.121794871794872,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0.121794871794872, 0, 0, 0, 1,
0, 1, 0, 1, 0.121794871794872, 0, 0, 0, 0.121794871794872, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.121794871794872,
1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0.121794871794872, 0, 0, 0,
1, 0, 0, 0, 0, 0, 1, 0, 0, 0.121794871794872, 0.121794871794872,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.115384615384615, 0, 0.115384615384615,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0.115384615384615,
0, 0, 0, 0.115384615384615, 0, 0, 0, 1, 0, 0, 0, 0.115384615384615,
0.115384615384615, 0.115384615384615, 0, 0, 0, 0, 0.115384615384615,
0, 0, 0.115384615384615, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0.115384615384615,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.115384615384615, 0, 0, 0, 0,
0, 0, 0, 0, 0.115384615384615, 0, 0, 0, 0.115384615384615, 0,
0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0.115384615384615,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.115384615384615, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.115384615384615, 0.115384615384615,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0.0506329113924051,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0506329113924051,
0, 0, 0, 0.0506329113924051, 0, 0, 0, 1, 0, 0, 0, 0.0506329113924051,
0.0506329113924051, 0.0506329113924051, 0, 0, 0, 0, 0.0506329113924051,
0, 0, 0.0506329113924051, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0.0506329113924051, 0, 0, 0, 0, 0, 0,
0, 0, 0.0506329113924051, 0, 0, 0, 0.0506329113924051, 0, 0,
0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.0506329113924051,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0506329113924051, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0.0506329113924051, 0.0506329113924051,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0.0961538461538462, 0, 0.0961538461538462,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.0961538461538462,
0, 0, 0, 0.0961538461538462, 0, 0, 0, 0, 0, 1, 1, 0.0961538461538462,
0.0961538461538462, 0.0961538461538462, 0, 0, 0, 0, 0.0961538461538462,
0, 0, 0.0961538461538462, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0961538461538462,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.0961538461538462, 0, 0, 0, 0,
0, 0, 0, 0, 0.0961538461538462, 0, 0, 0, 0.0961538461538462,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0961538461538462,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0961538461538462, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0961538461538462, 0.0961538461538462,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0,
1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0.0448717948717949, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0.0448717948717949, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0.0448717948717949, 0.0448717948717949, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0.0448717948717949, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.0448717948717949, 0, 0, 0, 0,
0, 0, 1, 0, 0.0448717948717949, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0.0448717948717949, 0.0448717948717949,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949, 0, 0.0448717948717949,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0.0448717948717949, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0.0448717948717949, 0.0448717948717949, 0, 1, 0, 0, 0.0448717948717949,
0, 0, 0.0448717948717949, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0.0448717948717949, 0, 0, 1, 0,
0, 0, 0, 0, 0.0448717948717949, 1, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0.0448717948717949, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(173L, 31L), .Dimnames = list(
    c("TCGA-AB-2803", "TCGA-AB-2805", "TCGA-AB-2806", "TCGA-AB-2807",
    "TCGA-AB-2808", "TCGA-AB-2810", "TCGA-AB-2811", "TCGA-AB-2812",
    "TCGA-AB-2813", "TCGA-AB-2814", "TCGA-AB-2815", "TCGA-AB-2816",
    "TCGA-AB-2817", "TCGA-AB-2818", "TCGA-AB-2819", "TCGA-AB-2820",
    "TCGA-AB-2821", "TCGA-AB-2822", "TCGA-AB-2823", "TCGA-AB-2824",
    "TCGA-AB-2825", "TCGA-AB-2826", "TCGA-AB-2828", "TCGA-AB-2830",
    "TCGA-AB-2832", "TCGA-AB-2833", "TCGA-AB-2834", "TCGA-AB-2835",
    "TCGA-AB-2836", "TCGA-AB-2837", "TCGA-AB-2838", "TCGA-AB-2839",
    "TCGA-AB-2840", "TCGA-AB-2841", "TCGA-AB-2842", "TCGA-AB-2843",
    "TCGA-AB-2844", "TCGA-AB-2845", "TCGA-AB-2846", "TCGA-AB-2847",
    "TCGA-AB-2848", "TCGA-AB-2849", "TCGA-AB-2851", "TCGA-AB-2853",
    "TCGA-AB-2854", "TCGA-AB-2855", "TCGA-AB-2856", "TCGA-AB-2857",
    "TCGA-AB-2858", "TCGA-AB-2859", "TCGA-AB-2860", "TCGA-AB-2861",
    "TCGA-AB-2862", "TCGA-AB-2863", "TCGA-AB-2865", "TCGA-AB-2866",
    "TCGA-AB-2867", "TCGA-AB-2868", "TCGA-AB-2869", "TCGA-AB-2870",
    "TCGA-AB-2871", "TCGA-AB-2872", "TCGA-AB-2873", "TCGA-AB-2874",
    "TCGA-AB-2875", "TCGA-AB-2877", "TCGA-AB-2879", "TCGA-AB-2880",
    "TCGA-AB-2881", "TCGA-AB-2882", "TCGA-AB-2884", "TCGA-AB-2885",
    "TCGA-AB-2886", "TCGA-AB-2887", "TCGA-AB-2888", "TCGA-AB-2889",
    "TCGA-AB-2890", "TCGA-AB-2891", "TCGA-AB-2895", "TCGA-AB-2896",
    "TCGA-AB-2897", "TCGA-AB-2898", "TCGA-AB-2899", "TCGA-AB-2900",
    "TCGA-AB-2901", "TCGA-AB-2903", "TCGA-AB-2904", "TCGA-AB-2908",
    "TCGA-AB-2909", "TCGA-AB-2910", "TCGA-AB-2911", "TCGA-AB-2912",
    "TCGA-AB-2913", "TCGA-AB-2914", "TCGA-AB-2915", "TCGA-AB-2916",
    "TCGA-AB-2917", "TCGA-AB-2918", "TCGA-AB-2919", "TCGA-AB-2920",
    "TCGA-AB-2921", "TCGA-AB-2924", "TCGA-AB-2925", "TCGA-AB-2927",
    "TCGA-AB-2928", "TCGA-AB-2929", "TCGA-AB-2930", "TCGA-AB-2931",
    "TCGA-AB-2932", "TCGA-AB-2933", "TCGA-AB-2934", "TCGA-AB-2935",
    "TCGA-AB-2936", "TCGA-AB-2937", "TCGA-AB-2938", "TCGA-AB-2939",
    "TCGA-AB-2940", "TCGA-AB-2941", "TCGA-AB-2942", "TCGA-AB-2943",
    "TCGA-AB-2944", "TCGA-AB-2946", "TCGA-AB-2948", "TCGA-AB-2949",
    "TCGA-AB-2950", "TCGA-AB-2952", "TCGA-AB-2954", "TCGA-AB-2955",
    "TCGA-AB-2956", "TCGA-AB-2959", "TCGA-AB-2963", "TCGA-AB-2964",
    "TCGA-AB-2965", "TCGA-AB-2966", "TCGA-AB-2967", "TCGA-AB-2969",
    "TCGA-AB-2970", "TCGA-AB-2971", "TCGA-AB-2972", "TCGA-AB-2973",
    "TCGA-AB-2975", "TCGA-AB-2976", "TCGA-AB-2977", "TCGA-AB-2978",
    "TCGA-AB-2979", "TCGA-AB-2980", "TCGA-AB-2981", "TCGA-AB-2982",
    "TCGA-AB-2983", "TCGA-AB-2984", "TCGA-AB-2985", "TCGA-AB-2986",
    "TCGA-AB-2987", "TCGA-AB-2988", "TCGA-AB-2990", "TCGA-AB-2991",
    "TCGA-AB-2992", "TCGA-AB-2993", "TCGA-AB-2994", "TCGA-AB-2995",
    "TCGA-AB-2996", "TCGA-AB-2998", "TCGA-AB-2999", "TCGA-AB-3000",
    "TCGA-AB-3001", "TCGA-AB-3002", "TCGA-AB-3005", "TCGA-AB-3006",
    "TCGA-AB-3007", "TCGA-AB-3008", "TCGA-AB-3009", "TCGA-AB-3011",
    "TCGA-AB-3012"), c("Offset", "CEBPA", "DNMT3A", "FAM5C",
    "FLT3_ITD", "FLT3_TKD", "IDH1", "IDH2", "KIT", "KRAS", "NPM1",
    "NRAS", "PHF6", "PTPN11", "RAD21", "RUNX1", "SMC1A", "SMC3",
    "STAG2", "TET2", "TP53", "U2AF1", "WT1", "complex", "minus5_5q",
    "minus7q", "plus8_8q", "plus21", "t_15_17", "t_8_21", "inv16_t16_16"
    )))

got an answer to start https://www.biostars.org/p/9533929/#9533942

A design matrix is a matrix of values of the grouping variable.
ANOVA needs such a matrix to know which samples belong to which group. Since limma performs an ANOVA it needs this design matrix.

You can create it using the model.matrix() method.
Afterwards you can tell limma which groups you want to compare.

For point 2 you define a contrast matrix defining the comparisons of interest by using the makeContrasts() method.
It is well explained in this manual's section 9.2


text:
Improving list speed: glycosylation example
I discussed a question with @gaspanic Python/Biopython - Replace amino acid residue on MSA with "z" from a list of unaligned positions
. The issue emerged was speeding up lists in Python. The solutions we discussed were

Raw code
Predefining a list size
deque

1. @gaspanic
my code
2. @gaspanic
my code
3. @gaspanic
my code

3.45-3.48
7.17-7.2
4.38-4.85
7.09-7.23
3.44-3.53
7.21-7.39

Run timings are seconds (3 replicated) above with @gaspanic's code being faster.
What improvements would speed this up?
Here I extended the data set
from itertools import repeat
def iterit(x):
    seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
    return ''.join(list(repeat(seq, x)))

if __name__ == '__main__':
    seqE3 = iterit(2800) # length is 201600
    seqE4 = iterit(28000) # length is 2016000
    seqE5 = iterit(280000) # length is 20160000

@gaspanic code
for seqbit in [seqE3,seqE4,seqE5]:
    seq = seqbit
    pos = [2, 5, 10, 11, 14]
    new_seq = []
    count = 0
    for s in seq:
        if s == "-":
            new_seq.append(s)
        else:
            count += 1
            if count in pos:
                new_seq.append('Z')
            else:
                new_seq.append(s)
   # print(''.join(new_seq))

My code:
# same variables and loop structure
newseq = list()
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        s = '-'
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
    newseq.append(s)
# print(''.join(newseq))

Predefining lists code (additional counter needed for @gaspanic's code:
# same variables and loop structure
new_seq = [] * len(seq)
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        s = '-'
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
    newseq[n] = s

Deque code simply any of the above substituting
newseq = list()

with
from collections import deque
newseq = deque()

Suggestions for efficiency improvements welcome.
A solution in Cython code would be cool. PyPI is not sought here.

I'd start by saying that Python is inherently slow and this feels like a case of perhaps unnecessary optimisation. With that in mind here are my suggested (untested) optimisations of @gaspanic code:
pos = {2, 5, 10, 11, 14} # make pos a set for O(1) member testing, and don't define inside the loop
for seq in (seqE3,seqE4,seqE5): # Define seq here rather than inside the loop, and use an immutable tuple instead of a mutable list
    new_seq = []
    count = 0
    for s in seq:
        if s == "-":
            new_seq.append(s)
        else:
            count += 1
            new_seq.append('Z' if count in pos else s) # use ternary operator

Thanks, the tuple did add speed, particularly tuple the iterator output. Using a set (the {} brackets) was a good solution and reduced runtime by ~0.5 seconds, which I'd not considered. It was better than tuple in this example. Repositioning pos adds a tiny amount, but agree its better practice and depends on the size of pos (in this case its small). Happy to discuss Python's slow speed in a summary.

Huge thanks to @Chris_Rands, this solution increases runtime by ~ 0.5 seconds uniformly across all coded solutions, so this code has a runtime of just over 3 seconds runtime on the benchmark machine. I incorporated my solution with @Chris_Rands solution and the final runtime I got was:
@gaspanic

1.7 - 1.8 seconds

This being twice the fastest speed in the question and four times the speed of my raw code.
My code

4.4 - 4.5 seconds

However a unique combination is needed:

tuple (Chris Rands)
set - important (Chris Rands)
an iterator
a list comprehension

In particular 3. and 4. are needed in combination. This is explained below.
def iterit(x):
    seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
    return tuple(''.join(list(repeat(seq, x))))

def mIterit(seqE3,seqE4,seqE5):
    pos = {2, 5, 10, 11, 14}
    for seqbit in (seqE3,seqE4,seqE5):
        seq = seqbit
        count = 0
        # loc = 0
        for s in seq:
            if s == "-":
                yield '-'
            else:
                count += 1
                if count in pos:
                    yield 'Z'
                else:
                    yield s

f __name__ == '__main__':

    seqE3 = iterit(2800) # length is 201600
    seqE4 = iterit(28000) # length is 2016000
    seqE5 = iterit(280000) # length is 20160000

    listComp = [i for i in mIterit(seqE3,seqE4,seqE5)]

The list comprehension (variable listComp) is the very last line and is crucial to the solution otherwise there is little increase in speed to @Chris_Rand's code. This alternative is to replace the last line with:
newseq = list()
for i in mIterit(seqE3,seqE4,seqE5):
    newseq.append(i)

A simple:
for i in mIterit(seqE3,seqE4,seqE5):
    pass

Is the same speed as the list comprehension, but its of no value because the output isn't saved to anything.
I was impressed by the interaction between the iterator and list comprehension. The other notable influence on speed is the set described by @Chris_Rand.
How does it work? tuples do increase speed by removing some functional of a list (immutability). Here the gain is small. Placing the pos variable outside the loop speaks for itself (0.1 seconds gain). set is interesting.
The iterator spews out a stream of results via yield which is where it attains its speed. The problem is capturing that stream outside the iterator is expensive. However list comprehensions are faster than standard loops and in this scenerio carry little/no detectable overhead. I was impressed. Up to this point I thought list comprehensions were a bit faster and a lot 'cooler' than loops. I hadn't understood there was a synergy with an iterator.
Is it all worth it? Python can easily achieve faster run times than Java using Cython but I will place the discussion in the chat forum next week.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Tool for cross referencing a newly found MOTIF to database of known MOTIFS
I need a tool or a function I can use in my code (R, or Python) that I can cross reference a MOTIF against known MOTIFS, a function that will take as input a MOTIF (a probability weight matrix, PWM which is a 4 by m matrix for an m-length motif) and spit out the closest -n- known motif(s). Is there such a function/tool or something close to it? Thanks much.

A popular tool for doing this is tomtom, part of the meme package.
You can try it via interactive web form here.
Once installed locally it can be run via command line
Additionally, it is available in the Bioconductor package memes as explained in Motif Comparison using TomTom

Thanks @malcook: TomTom seems like the tool.  The list of required libraries/tools to install first before tomtom can be installed is quite lengthy and even after that the installation is riddled with errors in the "make" step for me. Any hints on installation on Linux/Ubuntu?

Hi @ZebraFish - I have no hints - you are right: the meme suite is a bit of bear to install.
The R/BioConductor package 'universalmotif' also addresses this purpose and is  easier to install.  I have not (yet?) used it for such, but have been recently impressed by its documentation specifically that on "Motif comparisons and P-values".

Thanks @malcook: After much haggling and hacking around I discovered the simple way to install TomTom is actually through "conda install -c bioconda meme". This puts a slew of meme command line executables into a conda env. directory including TomTom. This simple method of installation was unnecessarily hard to find  if you are new to this field. Just looking at the meme suite webpage It's a curiosity why such simple tools (from a computational person's perspective motif search against a database of motifs is a very simple task) can be so obfuscated.

Answer from @m, converted from comment:
Hi @ZebraFish Motifs package Biopython
https://biopython.org/docs/1.76/api/Bio.motifs.html

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Simple template for calculating Shannon diversity index from FASTQ files with R
I have a bunch of files in FASTQ format and need to calculate the Shannon diversity index from all of them.
Is there a simple way to do that with R, some kind of step-by-step recipe?
My problem is not so much finding a function for calculating the index but how to do that on FASTQ files.
Grateful for every hint, Thank you!

A couple of questions: 1) Are you wanting to do the calculation over all the sequences in all the files together, or a separate calculation for each file? 2) Are your sequences already properly aligned, or does alignment need to be part of the solution? Please also note that in this context the Shannon index would only consider counts of different sequences, it would not consider how different the sequences are from each other (e.g. number of SNPs between two sequences doesn't matter, only that they're different).


text:
How do I read mmCIF files from Alphafold DB as coordinates to use in replicating their visualization for Unity VR?
https://alphafold.com/entry/Q9Y499
Starting at line 918 of the downloadable mmCIF file, it looks like there might be coordinates in x, y,z:
ATOM 1    N N   . TRP A 1 1   ? 12.549  -2.565  -21.811 1.0 96.06 ? 1   TRP A N   1 Q9Y499 UNP 1   W
ATOM 2    C CA  . TRP A 1 1   ? 12.477  -2.290  -20.365 1.0 96.06 ? 1   TRP A CA  1 Q9Y499 UNP 1   W
ATOM 3    C C   . TRP A 1 1   ? 11.221  -2.928  -19.812 1.0 96.06 ? 1   TRP A C   1 Q9Y499 UNP 1   W
ATOM 4    C CB  . TRP A 1 1   ? 12.469  -0.785  -20.103 1.0 96.06 ? 1   TRP A CB  1 Q9Y499 UNP 1   W
ATOM 5    O O   . TRP A 1 1   ? 10.206  -2.922  -20.501 1.0 96.06 ? 1   TRP A O   1 Q9Y499 UNP 1   W
ATOM 6    C CG  . TRP A 1 1   ? 13.795  -0.114  -20.270 1.0 96.06 ? 1   TRP A CG  1 Q9Y499 UNP 1   W
ATOM 7    C CD1 . TRP A 1 1   ? 14.286  0.508   -21.369 1.0 96.06 ? 1   TRP A CD1 1 Q9Y499 UNP 1   W
ATOM 8    C CD2 . TRP A 1 1   ? 14.791  0.098   -19.228 1.0 96.06 ? 1   TRP A CD2 1 Q9Y499 UNP 1   W

I am a bioinformatics noob (took an intro class in my undergrad CS degree at a Canadian university) but know some Unity game engine with which I want to build a standalone visualization tool for users with VR headsets.

"I want to build a standalone visualization tool" - Why, exactly? There are protein viewers in many languages (c - RasMol, Java - Jmol, C/Python - PyMol, Javascript = NGLViewer (?), etc) that will be able to show the structure in all sorts of detail.

more helpfully, I suppose the mmCIF format is here : http://ww1.iucr.org/iucr-top/cif/mm/index.html ; also yes, those look like the x/y/z coords (or Cartn_x, Cartn_y, Cartn_z) Again, there are mmCIF readers https://mmcif.wwpdb.org/docs/software-resources.html

@gilleain how about in VR?

Ok, that is slightly different, I see. For Jmol, at least, this was done a while ago - https://github.com/Jmol-OVR - in general, I would hope that some VR software could integrate with an external library. If you like, you could edit your question to focus on VR and I could make this comment into a proper answer

There are a few VR solution out there on the market, e.g. Nanome. The rate limiting step is the cost of the VR headset and user training so only biotech startups have them!

Don't
A protein is not only list of atomic coordinates. Bar for a proof of concept (i.e. learn to code in Unity VR), it is an endless spiral of hurt in a field where there is a lot of competition and were you to push a product you'd need to invest more manhours in publicity than in coding...
Connectivity
Here is the short of it. Chemistry is classically a graph network: nodes=atoms, edges=bonds. The connectivity of the atoms is dictated by their residue's 3-letter name, which is either

assumed —software know the 20 AAs (and several modifications, eg. SEP, 4x2 bases and some ligands, e.g. HOH— or
a special entry in the PDB or mmCIF file —in PDB it's a CONECT entry.
E.g. https://www.rcsb.org/ligand/CFF

AlphaFold2 does not have ligands, AlphaFill does. But scientists (not PR specialists) use PDB entries or run their own modelling. AlphaFold2 from EBI is rarely used: oligomers with ligands is where it's at. Especially if doing a drug discovery campaign (and red biotech can and does afford VR headsets).
Format
mmCIF is a deposition standard, PDB is the workhorse standard. mmTF is a cool idea that is not catching on (cf XKCD comic strip about standard proliferation...). Using PDB is easier: https://www.wwpdb.org/documentation/file-format-content/format33/v3.3.html
Representation
Problems arise when offering representations —everyone has their favourite. See NGL gallery for an array: https://nglviewer.org/ngl/gallery/index.html
Surface display has additional issues —PyMOL simply and openly reuses the Advance Poisson-Boltzmann Solvation tool as it's tricky.
Dragons in the minefield
In a PDB entry not from AF2, you can have disulfides (SSBOND and/or CIZ entry), isopeptide bonds (LINK), missing atoms, the UNK, UNX and UNL residues, alphatraces, insertion sequences, alt. occupancy, gaps, non-standard usages of all of these etc etc. and my favourite implied proximity bonding.

Okay, I get it. AlphaFold is not at all popular, heavily overhyped and has severe technical shortfalls. Good to know, thank you.

@M__ The algorithm (and its variants) are utterly-unutterably–undeniably phenomenal. It makes great complexes, thus making CoIP Westerns a thing of the stone age, does away with C-alpha traces in cryoEM and much more. It's the off-the-shelf EBI-AlphaFold2 models that are not as revolutionary as the publicity. The latest Trembl release is rather pointess in fact. Say one wants to see how a given pyrethrin-analogue drug binds in different arthropods? One would need to model the GABA channel as an oligomer. Norleucine insesitive Eco MetK? It's a missense in a patent not in Trembl. etc.

Thanks @MatteoFerla thats insightful. Pyrethrin-analogue in insects - this I do understand, E. coli metK dunno. CoIP Westerns being obsolete is a big statement - I did hear about the antibody binding prediction: I need to look into it in more detail and might forward this as a question in a month or so. Thanks again.

So basically "Don't bother for the above reasons" is what I'm hearing, is that correct?

Conditionally: if it's for learning, do go ahead but don't let feature crawl get you!
If it's to make a business that will compete against Nanome, yes, don't bother.

On the AF2 modelling side discussion, there is currently an explosion of variants —there's a repo https://github.com/sacdallago/folding_tools that's got things helpfully collected. However, ColabFold is the go to —although it does occasionally fail and demand MSA tinkering.


text:
Refactoring a script mapping ID to sequence

I wrote simple Python script which has an excessive run time, notably when using large data sets. For example a data set of 1 000 000 sequences.
I am seeking assistance refactoring the code and would be very grateful for any assistance.
from Bio import SeqIO
import re

def move_item(old, new):
    k, v = old.popitem()
    new[k] = v

def map_ids_to_aa_sequence_list(id_):
    id_ = id_.strip()
    for s in aa:
        if re.search(id_,s, re.I):
            return s

file = open(r'list_of_ids', 'r').readlines()
names = [for name.strip() in file if len(name) > 1]

handle = open('sequences.fasta','r')
sequences = list(SeqIO.parse(handle, 'fasta'))

seq_dict = {}

for s in sequences:
    seq_dict[str(s.description)] = str(s.seq)

aa = []

for k, v in seq_dict.items():
    sequence = '>' + k +'\n' + v + '\n'
    aa.append(sequence)

separated_aminoacid_sequences_set = set(map(map_ids_to_aa_sequence_list, names))

separated_aminoacid_sequences_list = list(separated_aminoacid_sequences_set)

with open('separated_sequences.fasta','w') as fs:
    for s in separated_aminoacid_sequences_list:
        fs.write(str(s))

On a brief glance the problem could be re but I'd need to look. Could you kindly look at this code here https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example
Are there any parallels?

Mapping millions of sequences to reference sequences is a common problem in bioinformatics, and has already been solved efficiently many times over. What is the problem that you are trying to solve? Are you creating your own code for educational purposes?

Hi @gringer. Thanks, the 'mapping' in the question I think (could be wrong) is the map() Python function. I edited the title - so lock of clarity is perhaps my error. The OP is using map() to try and speed up the function. (I personally didn't know whether yield would work better here). 'Recodifying' via map() is known practice. Lets see what the OP says. I'm happy to reedit accordingly.

The code is mapping amino acid sequences to reference sequences, but the ultimate goal of the problem is not clear from the question. Either way, more clarification is needed. If the goal is to make a hash lookup more efficient, then the question would be better answered in another StackExchange community. If the goal is to find some way to efficiently map amino acids to sequences, then plenty of other better existing solutions exist. Educational / learning is one reason to get answers about code optimisation.

This is beyond a single answer. 'Refactoring' (code efficiency, high performance code) is a huge topic. Like its own project.
The first thing I would personally absolutely target first is re here:
    if re.search(id_,s, re.I):

Python is not Perl and the re is EXTREMELY slow and in particular using search is especially slow, most notably when uncompiled (see notes and point 3 below).
How to avoid re.search
Could you revise this to a direct C style operation? This is the best solution if you get it working.
For example,

Using ==:

e.g.
if str(id_) == str(s):

if the case is important (like they're different) then...
if str(id_).lower() == str(s).lower():

You could dump the str() it ain't doing much its just bug insurance.
The other approach if you want to stick with re is to switch to match, that will be quicker.

Using match

e.g.
if re.match(id_, s, re.I):

match anchors  the search at the first character, i.e. id_[0], search goes everywhere. match is faster.

Sticking with search but combining with compile:

e.g.
id_ = id_.strip()
x = re.compile (id_)
if x.search(s, re.I): # don't use re. x replaces it

or
id_ = id_.strip()
x = re.compile (id_)
if x.match(s, re.I):

I have not put these in order of preference. My preference is for 1. but it will risk throwing logical errors (missing stuff).
If you are using re you must compile the search term. Its a pain but whilst compiled re is slow, uncompiled re is really, really slow. Disclaimer, I'm from a Perl background and from this perspective re is the pits because previously relied on the efficient regex (Perl).
The == and match could fall over if > is in one variable and not in the other, but the search will return the correct 'match' (i.e. not the output of 'match') at the cost of algorithm efficiency.
This isn't an exhaustive analysis of re but thats where I would start.
Please note this is a starting point from prior experience to pick a 'big issues'. It is useful for future, particularly for questions, to place time checks across the code to identify the bottlenecks.

Notes

The thing we don't know is what the fasta headers look like.
For every query in search there are $1.10^6$ fasta headers. So if you've $1.10^5$ queries thats $1.10^{11}$ searches for search thats why it can be taking up a lot of time.
Rewriting the code using Biopython objects throughout is good practice. SeqIO in particular are cool. AlignIO are algorithmically expensive however.


text:
Problems in conducting batch effect correction
I am now trying to conduct batch effect correction for my data matrix, whose samples (human) have different gender and come from different sequencing platforms. I have tried Combat, Combat-seq, and limma. However, the outcomes are bad. It seems that Combat-seq is based on Negative binomial regression models, Combat is based on Gaussian distribution, and limma is based on the linear model. I wonder if my data violates their assumption, so I check my data and get this outcome. It seems that my data follows Gamma distribution. I am wondering what I could do in this situation.

The data I am working on is DNAse-seq data, which is a matrix with 16 samples. As the samples come from different platforms, I would like to correct the batch effect. I am also considering gender information as there are male and female samples. To sum up, my questions are:

Do we need to consider the distribution of the data when we apply the batch correction (limma, combat, and combat-seq)?
If we need to consider the distribution, which method do I need to correct the data which follows gamma distribution?
Are there methods designed for read counts? (I know there are many methods designed for surat objects, but do we have others that are more general?)

Thank you very much for all kinds of help.

Dear @KeXu thanks for your question. Just a reminder if the answer below addresses your question, please do remember to "accept" it, and/or upvote (you have the rep to do this).

The ComBat-seq method requires raw counts which can be modelled as a negative binomial. The other methods are typically used with normalized counts transformed to log2 scale. Obviously for any method that tries to regress a factor out of a dataset you have to make sure that none of the experimental variables are confounded with the batch variables that you want to regress. More details can be given if you edit the post to say what the data are, how the experimental setup is and what the factors to regress are.

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.

Thank you very much for your reply. The data I am working on is DNAse-seq data, which is a matrix with 16 samples. As the samples come from different platforms, I would like to correct the batch effect.


text:
Polygenic Risk score calculation using vcf files
I have about 750 vcf files and I am looking to generate a Polygenic Risk Score for each individual. I am planning to use either PLINK or LDPhred to generate a Risk Score. I have phenotype information of these individuals.
Can anyone verify if the polygenic risk score is calculated using individual vcf files or should I look into combining the vcf files?

What kind of polygenic risk score? For what diseases? What genes? Using which of the many approaches?

For the tools you mentioned, you should convert to BED format. Other tools will accept VCF. It should work either way, combined or not, though you may have to adjust your pipeline depending on which you prefer. I'd imagine a combined file would be easier to work with.


text:
Can I use a tool for highly variable gene selection for selecting highly variable "bins" from an ATACSeq dataset?
I'm working with ATACSeq data from multiple tissue/cell types. The data is binned in 1 megabase bins. I'd like to identify the bins that are "highly variable" across the different tissue types. I can't seem to find a tool to do this, like there exists for identifying highly variables genes in a gene expression setting (e.g Scanpy's scanpy.pp.highly_variable_genes). I was wondering if I could use the gene selection tools for my data, and if there are any other suggestions?

I don't think is a good idea to apply the gene selection method on the ATAC-seq genomic regions because of the high sparsity of ATAC-seq data.
A common way is to select the most frequent regions (like the Signac does), or to discard the house-keeping regions (like the SnapATAC does).
And you can also use the episcanpy, a epigenomics version of scanpy.

Thanks for your answer. I guess in my case the cells are grouped per cell type (the data is aggregated per cell type), so sparsity shouldn't be an issue. Given this, is it still inappropriate to use gene selection methods?

If you think the sparsity is not an issue for your data, then scanpy.pp.highly_variable_genes may be applicable. Just be aware that the default thresholds for mean and disp are for single-cell RNA-seq data.


text:
Bar Graph of Expression Data from Seurat Object
I am looking for a way to present some single-cell RNA sequencing data on expression of a certain gene. My dataset has 2 variables, cell type and condition. I'm looking to create a grouped aligned scatter plot bar graph of the expression of the gene, where the y-axis is expression levels, while X axis is the different cell types, with each cell type having 2 bars (1 for healthy, 1 for diseased). My dataset has 3 healthy and 3 diseased samples, but all of the data is integrated into a Seurat object. To first create an aligned scatter plot bar graph, what I did was generate a DotPlot for the expression of gene X in each sample, split by cell-type. Then I realized that DotPlot() prints out a data table with values including avg.exp, pct.exp and avg.scaled.exp. I used the avg.exp metric to plot the data, but I was posting this question to see if there are other ways I can access and present expression data from a Seurat object. Is the way I've presented accurate? I've posted a sample bar graph of what I'm looking for, and the main question is how I can gather the expression levels across cell types AND diseased or healthy condition.
Thanks!

Regarding the other answer with the below code:
data.df <- object@meta.data[, c("celltypes", "condition")]
data.df["value"] <- object$RNA@data["geneX", ]
data.df <- data.df %>% mutate(group = paste(celltypes, condition, sep=":"))
data.df <- data.df %>% filter(value > 0)
data.df.2 <- data.df %>% group_by(group) %>% summarise(sd=sd(value),value=mean(value))
X <- str_split(data.df.2$'group', ":", 2, simplify=TRUE)

data.df.2$'celltypes' <- X[, 1]
data.df.2$'condition' <- X[, 2]

ggplot(data.df.2, aes(fill=condition, y=value, x=celltypes)) +
    geom_bar(color="black", position="dodge", stat="identity") +
    geom_errorbar(
        aes(ymin=value, ymax=value+sd), width=.2,
        position=position_dodge(.9)
    )

I find the plot matches more the results of Seurat differential expression if you DO NOT perform the step of removing the cells with 0 counts (data.df <- data.df %>% filter(value > 0)). To avoid having huge error bars you can use standard error of the mean instead of sd.

In single-cell data a bargraph with error bar is not usually shown since the gene expression does not have a normal distribution due to drop-outs. A violin-plot is commonly used. In your case, first set the active.ident to cell types
seuratObj <- SetIdent(seuratObj , value = "Cell.Types")

Then call the violin plot
VlnPlot(seuratObj , features = "geneX", split.by = "condition")

If you still specifically need bargraph, you can get the output of the violin plot as a dataframe and then draw the barplot.

Following the answer of @Ajay.
You can add + geom_boxplot() after the command to turn the violin plot to a box plot.
As for the barplot, you can do it with these:
# the object is your seurat object
data.df <- object@meta.data[, c("celltypes", "condition")]
data.df["value"] <- object$RNA@data["geneX", ]
data.df <- data.df %>% mutate(group = paste(celltypes, condition, sep=":"))

# remove cells that do not express the gene
data.df <- data.df %>% filter(value > 0)

data.df.2 <- data.df %>% group_by(group) %>% summarise(sd=sd(value),value=mean(value))
X <- str_split(data.df.2$'group', ":", 2, simplify=TRUE)
data.df.2$'celltypes' <- X[, 1]
data.df.2$'condition' <- X[, 2]

ggplot(data.df.2, aes(fill=condition, y=value, x=celltypes)) +
    geom_bar(color="black", position="dodge", stat="identity") +
    geom_errorbar(
        aes(ymin=value, ymax=value+sd), width=.2,
        position=position_dodge(.9)
    )

Dropout is very usual in single-cell data. If you really want to show the expression with barplot, I would suggest removing the zero values or using denoised data for better visualisation.


text:
Trinity de novo transcriptome assembly: samtools view failed to add PG line to the header
I have assembled the transcriptome a plant species using Trinity.
Here is the command:
Trinity --samples_file ./sample_file.txt --seqType fq --max_memory 80G --CPU 30 --output ./trinity_outdir

I tried to get transcript abundance using the following command:
/home/programs/build/trinityrnaseq/util/align_and_estimate_abundance.pl \
  --transcripts ./Trinity.fasta \
  --seqType fq \
  --samples_file ./sample_file.txt \
  --thread_count 6 \
  --gene_trans_map ./trinity_outdir.Trinity.fasta.gene_trans_map \
  --est_method RSEM \
  --aln_method bowtie2 \
  --prep_reference \
  --coordsort_bam \
  --output_dir /rsem_outdir/

And I keep getting the following error:
[E::sam_hrecs_update_hashes] Duplicate entry "TRINITY...." in sam header

samtools view: failed to add PG line to the header
I will appreciate your suggestions.


text:
Seurat: Error in FetchData.Seurat(object = object, vars = unique(x = expr.char[vars.use]), : None of the requested variables were found:
I was trying to do subsetting via:
pbmc = subset(pbmc, subset = nFeature > 1000 & nCount > 4000)

R returns:
Error in FetchData.Seurat(object = object, vars = unique(x = expr.char[vars.use]), :
None of the requested variables were found:

and then nothing, so it says cannot find requested variables and doesn't tell me what's missing...
I'm very confused. Thank you so much for your help!

I've come across the same error which seems to be a common issue when using the subset() function with Seurat objects. Sometimes, I encountered the same error and couldn't figure out the reason behind it, checking the variable's names was not a remedy in my case as I specified the exact same name for subsetting as to what I could see in my object metadata. However, I found a workaround to subset cells in a Seurat object using a different approach.
Instead of using subset(), you can try the following:
seurat_obj_subset <- seurat_obj[, <condition to be met>]

For example, if you want to subset a Seurat object called 'pbmc' based on conditions like having more than 1000 features and more than 4000 counts, you can use the following code:
pbmc = pbmc[, pbmc$nFeature > 1000 & pbmc$nCount > 4000]

From the steps in the PBMC_3k tutorial, it looks like you need to suffix the name with the assay as well, separated by an underscore (in this case, RNA):
pbmc <- subset(pbmc, subset = (nFeature_RNA > 1000) & (nCount_RNA > 4000))

"Error in FetchData.Seurat(object = object, vars = unique(x = expr.char[vars.use]), :
None of the requested variables were found:"

As the error is pointing out here, you should check your variable names that you have used here. It must match the variable names in the pbmc object that you want to subset.


text:
STAR vs Bowtie2
What is the fundamental difference between STAR and Bowtie(2). Specifically, what is the difference in their final output (regardless of run-time differences, speed, memory usage etc.). Both seem to work on the same kind of input *.fastq files without discriminating what they are associated with and generate an output version of the *.fastq files. Would my downstream analysis be incorrect if I use Bowtie2 where STAR was supposed to be used, and vice versa? Thank you.

What is the fundamental difference between STAR and Bowtie(2).

STAR for mapping spliced (i.e. with introns) short RNA-seq reads against a genome. Bowtie2 for mapping short reads without splicing.

Would my downstream analysis be incorrect if I use Bowtie2 where STAR was supposed to be used (?)

Yes. If you use bowtie2 to map spliced short reads against a genome, STAR-based downstream analysis is likely to be wrong.

... and vice versa?

If you map genomic reads with STAR, the result will be suboptimal. Use the right tool for the right task.

The seeds are clustered together by proximity to a selected set of ‘anchor’ seeds in the 2nd phase of the STAR algorithm,  but this •clustering/stitching/scoring process• is not present in Bowtie2.

mapping by performing a two-step process


text:
RNASeq NMF clustering using what kind of expression unit?
I would need some help regarding certain points in a RNASeq gene expression analysis (sorry if there are stupid questions, its my first project and i have no capable supervision). My workflow is the following:

TMM normalization

obtaining CPMs

filtering CPMs that are under 1 in half the samples

obtaining log2(CPM) - after I perform this step I get Inf values, should I treat those values as NA and remove them or should I use CPM + 1 (I know that adding +1 introduces certain bias)?

apply Median Absolute Deviation and use different thresholds to subset data

apply NMF clustering and find cluster specific genes and samples

use the sample clusters to obtain DEGs (using DESeq2) and proceed with GSEA

Now for the questions.. Firstly, I tried to look for the best practices in scientific papers and on forums but they vary in a huge manner.
In papers, such as https://pubmed.ncbi.nlm.nih.gov/26909576/ the authors use log2(CPM) for NMF clustering analysis, while looking online at other suggestions, they all suggest using variance stabilizing transformation or z-scores etc. Now, what I'm trying to do is find subtype specific genes and significant DE genes based on those subtypes.
From what I understand, VST transformation is ideal but why do all the papers that I've found only use TPM/CPM for subsequent analysis such as NMF clustering (other papers only mention gene expression or I missed that)?
I've tried both ways, with log2(CPM) and with VST and I've found that VST clusters contain more genes and tend to be of equal size (both genes and samples clusters) while log2(CPM) clusters returns a cluster 3 times bigger than the second (both for genes and samples). Both VST and log2(CPM) clusters were selected by the best cophenetic coef and also samples had silhouette widths above 0.6.
Now, what should be the best approach considering that I'm expected to write a decent paper with this dataset.
Thank you in advance!
PS: all the above are done in R


text:
bash script for interproscan REST
hey i have the problem with bash script. It should run interproscan REST
Im green in bash scripting. I found this script on the web:
#!/bin/bash
i=1
waitevery=30

mkdir -p out

for j in $(find `pwd` -type f -name "*.fa")
do
    echo "Iteration: $i; File: $j"
    filename=$(basename "$j")
    python3 iprscan5.py \
        --goterms \
        --pathways \
        --email=blabla@protonmail.com \
        --outfile=out/${filename} \
        --outformat=tsv \
        --quiet \
        $j & (( i++%waitevery==0 )) && wait
done

out:
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtr.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xts.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtt.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtu.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtv.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtw.fa

Edit:
I think I forgot to add --sequence, but something was still wrong after it.

Also, I'm assuming you're saving the script as submit.sh and then running it. How are you running it exactly?

@RamRS wait will wait until all background jobs have finished. (( i++%waitevery==0 )) will be true whenever the current value of i modulo waitevery is 0, so this is a really neat trick to make the script wait for all background processes to execute every N processes. Quite clever, that.

First of all never do for foo in $(find ...). That is very bad practice as it can't handle even slightly odd file names like those with a space. See bash pitfall #1.
Next, you don't need pwd, just run find -type f -name "*.fa" which will search in the current directory by default and unlike find `pwd`... , won't break if your path contains whitespace.
That's just general best practice, the real issue here is that you are probably running the script like this:
sh /path/to/script.sh

Right? And you are using a Debian-based system like Ubuntu, right? On such systems, /bin/sh is actually not bash but a very basic POSIX compatible shell called dash and dash doesn't understand the (( i++%waitevery==0 )) syntax.
Since you have a shebang line (#!/bin/bash), you can simply make your script executable (chmod a+x /path/to/script.sh) and then run the script directly like this:
/path/to/script.sh

Or, alternatively, run it explicitly with bash instead of sh:
bash /path/to/script.sh

Either should work and if your file names and paths are sane, the issues I mentioned above won't be relevant. But to keep things nice and clean, I would also change the script to (note that this assumes you are using Linux and GNU find):
#!/bin/bash
i=1
waitevery=30

mkdir -p out

find -type f -name "*.fa" -print0 |
  while IFS= read -r -d '' j; do
    echo "Iteration: $i; File: $j"
    filename=$(basename "$j")
    python3 iprscan5.py \
        --goterms \
        --pathways \
        --email="blabla@protonmail.com" \
        --outfile=out/"$filename" \
    --outformat=tsv \
    --quiet \
    "$j" & (( i++%waitevery==0 )) && wait
done

I asked them how they're running 'submit.shfor the reason you picked precisely - they're probably usingsh` and hitting a POSIX compatibility wall.

I'd guess that the specific error that's coming up is because there are spaces missing from the wait line. This seems to be a problem in the source script, i.e.:
$j & (( i++%waitevery==0 )) && wait

should probably be:
$j & (( i++ % waitevery == 0 )) && wait

but for me (following @terdon's reproduction advice via dash -c "i=9; ((++i%5==0)) && echo 'worked!'") that only changes the thing that is not found from i++%waitevery==0 to i++, suggesting that the post-increment operator is not supported in the shell. Post-increment is supported in bash, suggesting that this script is being run (as @RamRS commented) in an unexpected way.
Try running the script directly:
./interpro_batch_submit.sh

Or explicitly using bash:
bash ./interpro_batch_submit.sh

... but see @terdon's answer for a better approach to doing this.


text:
Compare illumina sequencing output with reference genome
I am new to whole genome sequencing. I want to compare illumina whole genome sequencing of mutants of a particular bacteria with a genome (not annotated) of the original strain to look for what the mutations are, how would I do about doing that?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.


text:
How can I find the maximum percent identity between two sets of (unaligned) sequences?
As the title states, I've got two sets of unaligned amino acid sequences (~25k sequences in one, ~3k in the other). I want to find the minimum distance between a sequence in the first set and a sequence in the second, expressed in terms of percent identity (so, maximum percent identity). So far I've only been able to find R packages that will compute percent identity for two already aligned sequences, which would then allow me to do a double loop over both sets of sequences to find the shortest distance between them. But this seems woefully inefficient - is there an easier/faster way? I also tried playing around with mmseqs2 but couldn't tell if there was a way to use it for anything other than clustering.
To clarify - there is a large amount of variation between the sequences in question, including in overall length (ranging from 148 residues to 872). They don't all belong to one family of homologs or anything. I'm expecting many of them to have quite low identity with each other. But the constraints of the problem I'm working on require me to put a number on that identity, or at least an upper bound.
I'm not constrained to using R, it's just what I'm most familiar with. Solutions with command line tools or python would also be fine - I'm sure there are better tools out there that I'm not aware of.

Wait ... are the 5' ends unequal lengths between target and query, in addition this is a 'no indel' solution sought? If so thats different and 'a few lines of code' is likely not the best solution.

When I've been in 'possibly' comparable situations I've just used local Blast and parsed out the top hit using Python. It may not be the most optimal package (actually its not), but its easy and works, whilst NCBI's Blast support desk is efficient. Each query in turn is tested via a loop against the whole local database.

I'm not very familiar with blast command line tools. What blast subcommand would I use for this?

I'll post it next week @empiromancer if no-one else does. Its an occasional thing, can't recall it on the spot. From memory easy, its parsing code for the output that's the main issue. I can take >6 months to respond, so one week is pretty quick :-)

Try mmseqs2. It is largely a faster version of blastp. I haven't used mmseqs2 before but its documentation seems simple enough.

I actually did play around with mmseqs2 but I couldn't figure out how to make it calculate percent identity. I could only figure out the clustering features.


text:
Converting from VCF to PLINK while maintaining the SNP names from the VCF file
Just as the title says, I am working in PLINK 1.9 trying to convert the following VCF into a PLINK file (binary or not either is fine), while maintaining the SNP names from the VCF file. I am on Windows, so it would be easier for me to utilize PLINK for this. Does anyone know how to achieve this? The link below is the VCF file I am using.
I tried using the --vcf command by itself, along with using the --vcf command in conjunction with the --recode command. I also tried using the --double-id command to see if that would affect anything, but it did not.
http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz

do you mean you've tried it already and it didn't preserve the SNP names? If so, please show which command you used.

Could you please explain why you need to convert a VCF file to PLINK format? I assume it's for another program, as VCF files can be read directly in PLINK.


text:
Transcriptome analysis
I am trying to assemble reads belonging to two different readlength. Is it a valid way since I am looking for common genes among the species I am assembling.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

My work focuses on identifying common genes amongst different species of insects. I am doing transcriptome analysis to find these genes. However, I am assembling reads from different readlength using Trinity. I am unable to understand how to proceed next.


text:
Difference between pileup and mpileup in VarScan and samtools
I'm going to call variants with VarScan from a pileup files created with samtools. I realized that there are in general two major possibilities to call variants, pileup as well as mpileup. The VarScan manual says that pileup is for single-sample, mpileup for multi-sample calling. Does that mean that one should use pileup to analyze clones, and mpileup if whole populations have been sequenced? What is the difference between both, pileup and mpileup? What would be the consequences to use the "wrong" pileup?
Moreover, I realized that samtools formerly also had the option pileup. However, pileup was later on deprecated. Thus, would you agree that using mpileup would do the job regardless of the sample?
Thanks a lot in advance!


text:
How to annotate list of proteins with KEGG orthology groups?
I've got a list of proteins, complete with Uniprot IDs (for some of them) and gene names. I'd like to group them into orthology groups based on data from KEGG, but I have no idea how to add a column with KEGG groups. I've tried poking around in the documentation for UniProt, but all I can find is this page on KEGG - no mention of how to link to KEGG data. I've also found these KEGG mapper tools (web portal and command line) but those don't seem to be what I need in order to cluster my proteins into groups of KEGG orthologues.
I've tried using R via the following:
library(KEGGREST)

kegg_ids <- KEGGREST::keggConv(
  target = "genes",
  source = paste("uniprot", uniprot_ids, sep = ':'))

But this turns up matches for only a fraction of my protein IDs... and anyway, only about half of them have those IDs - the rest are human proteins for which I only have the gene name.


text:
How are BLAST record scores calculated in biopython
Can someone explain how the score is calculated for Bio.Blast.Record.Description?  For example, if Record.Description.num_alignments is >1, is Record.Description.score simply the sum of the alignment scores for the description's hit?  Can alignments overlap?  How does the score of a given alignment deal with the scores of overlapping HSPs in that alignment?
I'm getting the Blast.Record object via NCBIXML.parse(), and would like to rely on the description data as much as possible to remove the need for code that parses the Record.Alignments objects, but I'm a nervous including things when I don't know exactly how they were created.
Cheers!
Tim

The results you'd get are the same as with any other NCBI BLAST+ tools, python does not compute it.
Biopython can run a blast query qblast on the internet or locally. To make it run locally you need to install blast tools on your machine and biopython will launch it as a subprocess. Here you can find the pertinent information: https://biopython.org/docs/dev/Tutorial/chapter_blast.html#running-blast-locally
Filtering data is as much a burden as it is important. Depending on your question you want to answer. Some people prefer a low e-value. The score is however affected by the length of the sequence since the number of matches is bigger, more details here: https://www.nlm.nih.gov/ncbi/workshops/2023-08_BLAST_evol/blast_score.html


text:
Read PDB file, extract dihedral angles, modify dihedral angles, reconstruct Cartesian coordinates, and write PDB file
As the title summarizes, I am trying to:

Read a PDB file (for example, 1enh.pdb).
Extract the backbone dihedral angles (phi, psi).
Modify the dihedral angles (phi, psi) arbitrarily (for example, replace them with np.random.uniform(-np.pi, np.pi, size = (56, 2)) in 1enh, since there are 56 pairs of dihedrals).
Reconstruct the Cartesian coordinates of the modified backbone.
Write a PDB file based on the modified backbone.

From this question it seems the tools for what I want are available in BioPython. There is also a set_angle() method. But I only have a very limited knowledge of Python and BioPython. I am stuck at Steps 3 and 4, and dubious on Step 5.
This is what I have tried so far:
#!/usr/bin/python

import numpy as np
from Bio.PDB import PDBParser, PICIO, PDBIO

## Step 1

parser = PDBParser()
structure = parser.get_structure("", "1enh.pdb")

## Step 2

structure.atom_to_internal_coordinates()
chain = list(structure.get_chains())[0]
ic_chain = chain.internal_coord
dihedrals = ic_chain.dihedra

## Step 3 - ??

# Modifiy (phi,psi). For example, replace (phi,psi) with
newdihedrals = np.random.uniform(-np.pi, np.pi, size = (56, 2))

## Step 4 - ??

modified = PICIO.read_PIC('internal')
modified.internal_to_atom_coordinates()

## Step 5

io = PDBIO()
io.set_structure(modified)
io.save('1enh_modified.pdb', preserve_atom_numbering = True)

Many thanks in advance.
EDIT: I do not care about possible clashes if the side-chain is included in the modified PDB. In any case, writing a modified PDB file that only includes the backbone is also fine.

POINT 5 READS : Write a PDB file based on the modified backbone. --> Just C-alpha or side-chains too ?? How to avoid clashes in case of the latter ?

Thanks for the comment. C-alpha would suffice. If side-chains are included, I do not really care about clashes, as I could determine them later and reject the modified structure.

Any Python package that checks PDB geometry  just for clashes ??

In PyMOL I remember it is possible to detect clashes, I do not know other functionality beyond it. In any case, my main interest is on tweaking the (phi,psi) angles and reconstruct the protein backbone in a PDB file. Side-chains could be present or not. Thanks.

ok, here my input pdb, example_short.pdb  :
ATOM      1  N   MET A   1      14.067  31.664  -1.639  1.00 80.87           N
ATOM      2  HN1 MET A   1      14.250  31.259  -0.699  1.00  0.00           H
ATOM      3  HN2 MET A   1      14.891  32.220  -1.945  1.00  0.00           H
ATOM      4  HN3 MET A   1      13.229  32.279  -1.596  1.00  0.00           H
ATOM      5  CA  MET A   1      13.824  30.564  -2.612  1.00 79.17           C
ATOM      6  HA  MET A   1      14.702  29.921  -2.671  1.00  0.00           H
ATOM      7  C   MET A   1      13.556  31.158  -3.989  1.00 78.52           C
ATOM      8  O   MET A   1      13.301  32.365  -4.085  1.00 79.40           O
ATOM      9  CB  MET A   1      12.634  29.724  -2.148  1.00 83.38           C
ATOM     10  HB1 MET A   1      12.600  29.755  -1.059  1.00  0.00           H
ATOM     11  HB2 MET A   1      11.725  30.170  -2.552  1.00  0.00           H
ATOM     12  CG  MET A   1      12.680  28.264  -2.580  1.00 85.71           C
ATOM     13  HG1 MET A   1      12.066  28.145  -3.473  1.00  0.00           H
ATOM     14  HG2 MET A   1      12.274  27.650  -1.776  1.00  0.00           H
ATOM     15  SD  MET A   1      14.354  27.697  -2.951  1.00100.24           S
ATOM     16  CE  MET A   1      15.137  27.675  -1.338  1.00 97.10           C
ATOM     17  HE1 MET A   1      15.118  28.678  -0.913  1.00  0.00           H
ATOM     18  HE2 MET A   1      14.597  26.991  -0.683  1.00  0.00           H
ATOM     19  HE3 MET A   1      16.170  27.342  -1.441  1.00  0.00           H
ATOM     20  N   GLU A   2      13.604  30.306  -5.018  1.00 75.32           N
ATOM     21  H   GLU A   2      13.628  29.291  -4.794  1.00  0.00           H
ATOM     22  CA  GLU A   2      13.626  30.684  -6.444  1.00 72.66           C
ATOM     23  HA  GLU A   2      13.508  29.732  -6.961  1.00  0.00           H
ATOM     24  C   GLU A   2      14.903  31.384  -6.843  1.00 65.44           C
ATOM     25  O   GLU A   2      15.489  31.127  -7.892  1.00 67.46           O
ATOM     26  CB  GLU A   2      12.516  31.672  -6.810  1.00 72.71           C
ATOM     27  HB1 GLU A   2      12.094  32.068  -5.886  1.00  0.00           H
ATOM     28  HB2 GLU A   2      12.960  32.486  -7.383  1.00  0.00           H
ATOM     29  CG  GLU A   2      11.397  31.081  -7.622  1.00 77.72           C
ATOM     30  HG1 GLU A   2      10.901  31.877  -8.178  1.00  0.00           H
ATOM     31  HG2 GLU A   2      11.809  30.354  -8.321  1.00  0.00           H
ATOM     32  CD  GLU A   2      10.385  30.393  -6.731  1.00 79.89           C
ATOM     33  OE1 GLU A   2       9.171  30.574  -6.987  1.00 82.69           O
ATOM     34  OE2 GLU A   2      10.824  29.702  -5.779  1.00 72.73           O
ATOM     35  N   ASN A   3      15.194  32.399  -6.043  1.00 60.12           N
ATOM     36  H   ASN A   3      14.712  32.444  -5.123  1.00  0.00           H
ATOM     37  CA  ASN A   3      16.140  33.452  -6.363  1.00 56.21           C
ATOM     38  HA  ASN A   3      16.072  33.715  -7.419  1.00  0.00           H
ATOM     39  C   ASN A   3      17.548  32.968  -6.047  1.00 50.40           C
ATOM     40  O   ASN A   3      18.513  33.666  -6.331  1.00 50.95           O
ATOM     41  CB  ASN A   3      15.821  34.667  -5.493  1.00 57.10           C
ATOM     42  HB1 ASN A   3      15.913  34.381  -4.445  1.00  0.00           H
ATOM     43  HB2 ASN A   3      16.538  35.457  -5.717  1.00  0.00           H
ATOM     44  CG  ASN A   3      14.425  35.196  -5.732  1.00 62.07           C
ATOM     45  OD1 ASN A   3      14.220  36.001  -6.641  1.00 69.12           O
ATOM     46  ND2 ASN A   3      13.464  34.748  -4.923  1.00 57.21           N
ATOM     47 1HD2 ASN A   3      12.486  35.079  -5.045  1.00  0.00           H
ATOM     48 2HD2 ASN A   3      13.693  34.068  -4.170  1.00  0.00           H
ATOM     49  N   PHE A   4      17.643  31.764  -5.492  1.00 47.14           N
ATOM     50  H   PHE A   4      16.781  31.193  -5.378  1.00  0.00           H
ATOM     51  CA  PHE A   4      18.914  31.219  -5.038  1.00 46.77           C
ATOM     52  HA  PHE A   4      19.696  31.915  -5.340  1.00  0.00           H
ATOM     53  C   PHE A   4      19.239  29.888  -5.695  1.00 48.73           C
ATOM     54  O   PHE A   4      18.490  28.919  -5.594  1.00 49.50           O
ATOM     55  CB  PHE A   4      18.922  31.095  -3.507  1.00 42.32           C
ATOM     56  HB1 PHE A   4      18.091  30.456  -3.207  1.00  0.00           H
ATOM     57  HB2 PHE A   4      19.862  30.635  -3.201  1.00  0.00           H
ATOM     58  CG  PHE A   4      18.787  32.415  -2.806  1.00 39.42           C
ATOM     59  CD1 PHE A   4      17.545  32.878  -2.408  1.00 40.76           C
ATOM     60  HD1 PHE A   4      16.673  32.232  -2.509  1.00  0.00           H
ATOM     61  CD2 PHE A   4      19.881  33.258  -2.648  1.00 35.13           C
ATOM     62  HD2 PHE A   4      20.857  32.941  -3.015  1.00  0.00           H
ATOM     63  CE1 PHE A   4      17.397  34.138  -1.889  1.00 36.49           C
ATOM     64  HE1 PHE A   4      16.400  34.520  -1.669  1.00  0.00           H
ATOM     65  CE2 PHE A   4      19.757  34.474  -2.043  1.00 32.88           C
ATOM     66  HE2 PHE A   4      20.640  35.090  -1.873  1.00  0.00           H
ATOM     67  CZ  PHE A   4      18.507  34.929  -1.643  1.00 43.95           C
ATOM     68  HZ  PHE A   4      18.401  35.892  -1.144  1.00  0.00           H
END

here my code :
from Bio import PDB

from Bio.PDB import PICIO, PDBIO

from typing import TypedDict, Dict, Tuple

parser = PDB.PDBParser(PERMISSIVE=1, QUIET=1)

structure: PDB.Structure.Structure = parser.get_structure("prova", "example_short.pdb")

structure.atom_to_internal_coordinates()
chain: PDB.Chain.Chain = list(structure.get_chains())[0]
ic_chain: PDB.internal_coords.IC_Chain = chain.internal_coord
d:  Dict[Tuple[PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey],
                 PDB.internal_coords.Dihedron] = ic_chain.dihedra

cnt = 1
for key in d:

    if key[0].akl[3] == 'N':
        if key[1].akl[3] == 'CA':
            if key[2].akl[3] == 'C':
                if key[3].akl[3] == 'N':

                    print ('\n',cnt,' :   ',  [x.akl[3] for x in key], d[key].angle)

                    d[key].angle += 45

                    cnt += 1

cnt = 1
for key in d:

    if key[0].akl[3] == 'N':
        if key[1].akl[3] == 'CA':
            if key[2].akl[3] == 'C':
                if key[3].akl[3] == 'N':

                    print ('\n',cnt,' :   ',  [x.akl[3] for x in key], d[key].angle)

                    cnt += 1

structure.internal_to_atom_coordinates(verbose = True)

io = PDBIO()

io.set_structure(structure)

io.save('atom_coord.pdb',  preserve_atom_numbering=True)

here an image of the original vs new quadripeptid (cartoon repr, just backbone) :
original is rainbowed one (blue to red), magenta is the new one
saved as atom_coord.pdb

I believe I changed φ (Phi) adding 45° since

Now, a φ angle is N-CA-C-N

It is worth mentioning that changing
 `d[key].angle += 45` to `d[key].angle += 6645`

works too. Apparently, need confirmation from someone more knowledgeable,
`structure.internal_to_atom_coordinates`

can cope with more than 360° angle values.

A nice and simple experiment is making a spaghetti (spaghetto?) into a helix with φ=-57.8 and ψ=-47.0

For the record, I wanted to award the bounty to this answer, but it was auto-assigned :(

I do not know if there is a function in Biopython.PDB for dihedral space transformations, but due to the lack of answers I would say no. Biopython.PDB is nice and clean, but I would use something with more molecular mechanics power... or should I say potential (:drum::laughing:).
Personally, I use PyRosetta, but it is not open source (code available to academics), so others may opt for OpenForceField.
PyRosetta requires a licence (free for academia), which is different from the Rosetta proper licence.
The code below uses a helper module I wrote, so I'll use the installer therein (which is shorthand for the official pip wheel, which provides different options etc., while the official shorthand is actually pyrosettacolabsetup)
pip3 install pyrosetta-help
install_pyrosetta -u  -p

Ligands?
First lets make sure there's nothing interesting ligand wise.
NB. RCSB has a new shiny API which may be better... Not tried it: currently using European.
import requests

def retrieve_nonpolymer_data(entry_id: str) -> dict:
    """cf. https://www.ebi.ac.uk/pdbe/api/doc/"""
    response: requests.Response = requests.get(f' https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/{entry_id.lower()}')
    response.raise_for_status()
    return response.json()[entry_id.lower()]

info = retrieve_nonpolymer_data('1enh')
chem_resns = [c for entity in info for c in entity.get('chem_comp_ids', [])]
print(chem_resns)  # ['HOH']

If is were something that is not water, I made a list here, which is based off Peter Curran's repo for the HotSpots algorithm, which allows the elimination of worthless ligands like glycol.
I mention this because the route differs if there is a ligand.
Init and read
import pyrosetta
import pyrosetta_help as ph
from types import ModuleType

logger = ph.configure_logger()
pyrosetta.init(extra_options=ph.make_option_string(no_optH=False,
                                                ex1=None,
                                                ex2=None,
                                                #mute='all',
                                                ignore_unrecognized_res=True,
                                                load_PDB_components=False,
                                                ignore_waters=True)
                               )

# shorthand for sanity
prc: ModuleType = pyrosetta.rosetta.core

# this route strips ligands
pose: Pyrosetta.Pose = pyrosetta.toolbox.pose_from_rcsb('1enh')

Note ignore_unrecognized_res, load_PDB_components and ignore_waters.
Angles
This will not happen to an non-Rosetta user, but I get tripped out by the fact that ψ & φ data is not stored in residue, but in the pose.
for residue in pose.residues:  #: prc.conformation.Residue
    # residue.get_phi() # ---> Error!
    pass

For the residue index in pose, Rosetta was written in Fortran so is one-indexed.
import numpy as np
import numpy.typing as npt
phi:npt.NDArray[float] = np.fromiter(map(pose.phi, range(1, pose.total_residue()+1)), dtype=float)
psi:npt.NDArray[float] = np.fromiter(map(pose.psi, range(1, pose.total_residue()+1)), dtype=float)

Change φ/ψ angles
Time for some fun in dihedral space, say:
r = 2
phi = 180
pose.set_phi(r, phi)

For modelling, I wrote some code for setting spans with the archetypical angles for different SS: here. I do not recall whence the data came (maybe Wikipedia?), but it may be handy anyway.
Save
To save a PDB file:
pose.dump_pdb('.pdb')

There is also pose.dump_scored_pdb, but that is only for friends of J. Willard Gibbs.
Minimisation
It is said in the comments that energy minimisation is not required. But here is a snippet that minimises the sidechains only
scorefxn = pyrosetta.get_fa_scorefxn()  # ref2015 will do.
before: float = scorefxn(pose)
cycles = 5
movemap = pyrosetta.MoveMap()
movemap.set_bb(False)  # backbone
movemap.set_chi(True)  # sidechain
relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxn, cycles)
relax.set_movemap(movemap)
relax.apply(pose)
after: float = scorefxn(pose)
print(after - before)  # Willard says "Remember free energy is a potential so it should be negative"

If the BB tweaks are minimal, there's also the option of scoring/minimising against the electron density etc.

Thanks for the answer. In BioPython there is functionality for setting angles in the internal_to_atom_coordinates() method. To test your solution, could you edit to a single script that does Steps 1-5 and sets all the phi/psi angles? I am not familiar with PyRosetta. Something that worries me: are you doing some behind-the-scenes optimization? I just want to make the transformation of going from dihedral to Cartesian space. A test of correctness: get back the modified dihedrals from the reconstructed atom coordinates.

Uhm... It does strike me as odd that you cannot alter the dihedrals in that case —let me check. Bio.PDB is rather basic so I would assume it might be convoluted hence the lack of Bio.PDB replies.
PyRosetta does add protons during reading and strips waters if told to do so. If told to autoload the chemical components (ligands) and the structure has an anomaly in its ligand (i.e. it is not what it ought to be) it will not be distorted on load unless minimised.

I answered about doing it with PyRosetta as I assumed that as BioPython PDB module is basic and the fact it was being asked made me think it was not possible.
Turns out I was wrong. So I am making a separate answer —as I do not actually believe doing it with Bio.PDB is the best solution.
The provided code from above with type hints is:
import numpy as np
from Bio import PDB

from typing import TypedDict, Dict, Tuple

parser = PDB.PDBParser()
structure: PDB.Structure.Structure = parser.get_structure("", "/Users/matteo/Downloads/1enh.pdb")
structure.atom_to_internal_coordinates()
chain: PDB.Chain.Chain = list(structure.get_chains())[0]
ic_chain: PDB.internal_coords.IC_Chain = chain.internal_coord
dihedrals:  Dict[Tuple[PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey],
                 PDB.internal_coords.Dihedron] = ic_chain.dihedra

I got the types simply by using type and the attributes of note with dir. The latter helps me explore the Bio.PDB.internal_coords.Dihedron telling me of the angle attribute:
d: PDB.internal_coords.Dihedron = list(dihedrals.values())[0]
d.angle: np.float64

A change here appears it the main PDB.Structure.Structure instance, i.e. it is not a disconnected copy.
The class PDB.internal_coords.AtomKey appears to be hard to instantiate by humans and may not work as expected —help is gibberish so it is a case for inspect.getsource. So iterating across the dihedrals for the Ramachandran dihedrals is the solution. Now, a φ angle is N-CA-C-N and a ψ is C-N-CA-C (in PDB atom names the carboxylic carbon is C ). When the __str__ method of one of these PDB.internal_coords.AtomKey instances is called some values are given, which are available in the tuple attribute akl
atom_key: PDB.internal_coords.AtomKey = list(ic_chain.dihedra.keys())[0][0]
print(atom_key)  #3_R_N
print(atom_key.akl) #('3', None, 'R', 'N', None, None)

The first entry is the residue number and is a string because of insertion codes.
The fourth entry is the atom name. So iterating across the keys for the correct sequence is the trick —I know there's bounty asking for a complete solution, but it's simple, so anyone care to fill out the blanks and get the meaningless points, please do!

Thanks. Indeed I am actually looking for a complete solution that takes care also of the 'meaningless points' and that has a test of correctness.

I am at a loss here, one question so biopython doesnt give you just the backbone/Ramachandran dihedrals but you need to sort them out ?


text:
Sequence format: grouping sequences within a flat file
Is there a sequence format allowing multiple sequences to be grouped  across a flat file? The specific application is stacking separate nucleotide sequences against a given RNA secondary structure.
To rephrase the same question I'm looking to stack viral nucleotide data in different groups across a single flat file. FastaQ provides an example:
Traditional fastQ,
@sequence_name
AGCTGAT
+
...(quality score)

@sequence2
AGCTGAC
+
...(quality score)

What I'm looking for is a format,
@Gp1_virus1-2-3
AGACGTAG
+
AGAGACAA
+
AGAGACAT

@Gp2_virus1-2
AGATGTAG
+
AGATACAA

In the above example there are two groups comprising 3 (group 1) and 2 (group 2)  individual sequences (viruses). Here 'group' represents different structures.
Stockholme format has potential. My understanding of this format, which is basic, is it stacks a sequence, together with structural information. What is sought is stacking multiple sequences and then a structure format. Any clarification would help.
There are multiple XML formats, which may be relevant.
Just to mention a 'flat file' is needed to allow manual inspection. Dumping JSON is a solution, but isn't friendly for optical inspection.

Application
The immediate downstream application is structure/sequence viewers. I think at least one. viewer will accept Stockholme 1.0 format. rFAM is very similar (@terdon's comment)
Summarised Question The question could be rephrased as to whether multiple RNA families can be stacked inside one Stockholme 1.0 or else rFAM format?
Some of the members here are likely familiar with pFAM. So could multiple protein structures be stored and retrieved inside one pFAM file.
@gringer's contribution is very welcome.

What is this for? I mean, since you don't seem to need something that will be understood by a specific downstream program, can't you just define your own format?

Thanks, @terdon I've updated the post. A hybrid solution may be needed (last resort).

An unmapped SAM file format allows for read groups. You can use samtools import to create unmapped SAM files from FASTQ files:
samtools import -R virus-1-2 in12.fq > 12_reads.sam
samtools import -R virus-1-2-3 in123.fq > 123_reads.sam
samtools merge out.sam 12_reads.sam 123_reads.sam

Perfect, thank you. It can be manipulated via Biopython SamtoolsCalmdCommandline

I'm awarding 'solved' to @gringer, because thats the context of the question.
The best solution I could find is stacking Stockholm (sth)  formats together. In this example the RNA families are associated with one structure:
# STOCKHOLM 1.0
#=GF ID    UPSK
AF035635.1/619-641             UGAGUUCUCGAUCUCUAAAAUCG
M24804.1/82-104                UGAGUUCUCUAUCUCUAAAAUCG
J04373.1/6212-6234             UAAGUUCUCGAUCUUUAAAAUCG
M24803.1/1-23                  UAAGUUCUCGAUCUCUAAAAUCG
#=GC SS_cons                   .AAA....<<<<aaa....>>>>
//

Biopython allows sth format to be input.
from Bio import AlignIO
from pathlib import Path

file = 'test.sth'
path = '/Volumes/mypath/testSth'
alignment = AlignIO.read(open(Path(path, file)), "stockholm")
print("Alignment length %i" % alignment.get_alignment_length())
print (alignment.column_annotations['secondary_structure'])

Output,
Alignment length 23
.AAA....<<<<aaa....>>>> # this is the RNA structure

If we look the alignment object via
from pprint import pprint
pprint(vars(alignment)) # or via print(alignment.__dict__)

The output

{'_annotations': {},  '_per_col_annotations': {'secondary_structure':
'.AAA....<<<<aaa....>>>>'},
'_records':
[SeqRecord(seq=Seq('UGAGUUCUCGAUCUCUAAAAUCG'),
id='AF035635.1/619-641', name='AF035635.1',
description='AF035635.1/619-641', dbxrefs=[]),
SeqRecord(seq=Seq('UGAGUUCUCUAUCUCUAAAAUCG'), id='M24804.1/82-104', name='M24804.1', description='M24804.1/82-104',
dbxrefs=[]),
SeqRecord(seq=Seq('UAAGUUCUCGAUCUUUAAAAUCG'), id='J04373.1/6212-6234', name='J04373.1',
description='J04373.1/6212-6234', dbxrefs=[]),
SeqRecord(seq=Seq('UAAGUUCUCGAUCUCUAAAAUCG'), id='M24803.1/1-23', name='M24803.1', description='M24803.1/1-23',
dbxrefs=[])],  'annotations': {}} {'secondary_structure':
'.AAA....<<<<aaa....>>>>'}

The reason vars(alignment) is required is to internally construct an alignment object in Stockholm format, here's the first part of the skeleton 'pseudo-code',
from Bio.Align import MultipleSeqAlignment
from Bio.SeqRecord import SeqRecord

outAlign = MultipleSeqAlignment([])
record in alignment:
    name = record.id
    myseq = record.seq
    sr = SeqRecord(Seq(name), myseq, '', '')
    outAlign.append(sr)

The above code constructs a basic alignment, the bit thats missing is to add the following key/value to the alignment object of {'secondary_structure': '.AAA....<<<<aaa....>>>>'}.
Its complicated, but you get use to it.

https://github.com/biopython/biopython/blob/master/Bio/AlignIO/StockholmIO.py


text:
What is the most accurate approach for de Novo sequencing?
I'm trying to decide between PacBio HiFi or Illumina sequencing platforms for sequencing the genome of aChrysina scarab. We want to identify the color pattern loci and also perform a phylogenetic analysis as we'll be sequencing several species of that genus.
When considering a hybrid approach I was told by someone at a sequencing facility that adding short reads to polish the HiFi reads assembly lowers the quality of the assembly and introduces more errors that corrections. Why is this? I haven't been able to find this statement in the literature.

HiFi reads can resolve most repeats in a genome. Illumina reads from these repetitive regions are likely to get mismapped and these mapping errors may become consensus errors. Most short-read polishing tools developed in the old days (e.g. pilon) will reduce the contig consensus quality. Furthermore, when you have a non-inbred diploid genome, HiFi reads can locally phase a large portion of the genome, while Illumina reads are likely to cancel/discard most HiFi phasing information, which would lead to a worse and less informative assembly.
In theory, if you use Illumina reads the right way, you can further improve HiFi contig consensus like what the telomere-to-telomere team is doing here. However, this pipeline is for haploid/inbred genomes only. There are no tools for generic genomes. Developing your own is beyond the capability of most people in the field. To this end, the person at your sequencing facility is mostly right: don't polish HiFi assemblies with Illumina reads.

could you please expand on this: "reads are likely to smear the phasing"?, I'm not sure I got your point.

@Caternina See the revised answer.

Very interesting. Personally, I think the most important thing than anything I can input is it is worth asking about the HiFi quality index being considered and whether this was compatible with your biological objective.
Just to highlight that insect genomes are very different to genomics here - traditionally and even now - there has been a disproportionate focus on dipteran genomes. Okay certain Lepidoptera, namely the silk moth (Bombyx), is an advanced area of study. Beetles genomes are even more different still, agreed they holometabolous insects but still along way from dipterans. In addition, its worth pointing out beetles were late on the genome sequencing scene.
I don't know the precise biological objective - but I'm going to guess it.  I'm assuming your wanting SNP variation to associate with phenotype colour variation, for the target genes (just guessing).

HiFi reads

With regarding mixing the genome sequencing approaches: what is the quality scale being used? If its lowering the chromosomal assembly (traditionally a key difficulty in insects), or weakening contig assembly then ok I could imagine. This is very different DNA. If the quality index is relating to lowering SNP accuracy I would be surprised, for others to comment, but just can't see it.
I am aware that a key issue is that insect genome, specifically with reference to dipteran genomes, has caused extreme difficult assembling due to the high AT content and repeat regions, with chromosomal identification a traditiona nightmare and requiring wet lab chromosomal mapping (Diptera). I think beetles are an outlier even on this spectrum - against other insects. The advice you might have been given could be with this issue in mind, particularly if it was from another entomologist.
I don't know if this genome has been sequenced, the implication of the question (title) is that is hasn't been sequenced yet and which leads to the second issue,
JBS Haldane

God has an inordinate fondness for stars and beetles.

There are an awful lot of beetles and a lot of those genomes not sequenced yet.
Personally I think you don't want accurate (really to maximise) genome assembly. In that scenario its got to be Illumina short-read, because this is a SNP study, where high accuracy calling from read depth is needed, given my guessed objective.
If the genome hasn't been sequenced there could be pressure to enhance genome assembly rather than target your biological objective. This could be why there has been mixed advice, between what's good for the community and what's needed for your biological objective, this being colour variation genes - for very pretty beetles.
I think, personally, its not necessarily a technical issue its either a win for others, or win for you. The only caveat is if the genes of interest are some extensive, physically linked gene family, then of course short-read may not be cool. I am not aware of that is a key feature of insect genomes.

In insect genomes there has been misplaced excitement ... new insect genome ... hang on thats just unordered  contigs. This will happen using a short-read only approach.