# Question
Title: Microbial Classification with 16S Nanopore Data
I have recently been tasked with analyzing a heterogeneous microbial sample composed of 16S data generated on a MinION sequencer, with the goal of assigning taxonomies and calculating respective abundance.

It is my understanding that traditional preprocessing methods such as Dada2, Unoise3, and Deblur would likely not be appropriate for this kind of data given the poorer read quality inherent with nanopore sequencing. To give an example, the average read quality of this current experiment is a Phred score of 15 and an average read length of 1500bp.

Furthermore, it would seem best to try and utilize a classifier that makes use of a full length 16S reference, rather than just the V3 region, to try and offset the higher error rate. But are there any classifiers that do this out of the box?

I am very interested to hear of any advice or input the community has on this. Thanks!

# Answer
Use kraken2 with one of the 16S databases. Kraken2 has worked well for me with nanopore sequences. Pavian can be used to inspect the Kraken2 report files visually.

Using uncorrected called reads has worked for me. Any errors from nanopore sequencing should be inconsistent enough (i.e. not always lead towards the same wrong organism) that the noise doesn't matter.

Note: if you're using the MiniKraken database instead of a comprehensive 16S database, there might be some false positive species assignments depending on which kmers get pruned out.

Another possibility is centrifuge, but I'm not sure if that works for 16S amplicon sequencing; it's primarily for metagenomic analysis.
> 3 votes

# Answer
I would second the answer gringer gave recommending Kraken2. I would suggest using bracken to process the results from Kraken2 if you want to estimate species abundance (the raw Kraken2 results should not be used for this). I would recommend Pavian for visualising the results, it produces beautiful Sankey plots in a fairly pain-free manner. If you want to save yourself some work the 16S tool provided by Oxford Nanopore as part of epi2me produces good results for limited effort - although this is a bit of a proprietary black box.
> 2 votes

---
Tags: nanopore, minion, microbiology
---


# Question
Title: Clustering information saved in Seurat object
Seurat V2 had a option to find clustering information saved in object: PrintFindClustersParams(object = pbmc). How can I get the same clustering parameters from objects in Seurat3? Do I need to manually save clustering parameters or they are saved automatically?

Update: the command line `object@commands` works perfectly for `...$FindClusters`. However, if I would like to find information about dimensions used there is an error:

```
...$FindNeighbors.SCT.pca Error in [[.Seurat(x, i, drop = TRUE) : Cannot find 'FindNeighbors.SCT.pca' in this Seurat object

```

# Answer
I am not sure of the output of `PrintFindClustersParams`, as I haven't use it before, but try with `object@commands`. It returns a list containing every command and parameters used.

For example

```
...
$FindNeighbors.SCT.pca
Command: FindNeighbors(seurObj, dims = 1:20, verbose = vrb)
Time: 2019-12-11 12:41:12
reduction : pca
dims : 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
assay : SCT
k.param : 20
compute.SNN : TRUE
prune.SNN : 0.06666667
nn.method : rann
annoy.metric : euclidean
nn.eps : 0
verbose : FALSE
force.recalc : FALSE
do.plot : FALSE
graph.name : SCT_nn SCT_snn

$FindClusters
Command: FindClusters(seurObj, verbose = vrb)
Time: 2019-12-11 12:41:17
graph.name : SCT_snn
modularity.fxn : 1
resolution : 0.8
algorithm : 1
n.start : 10
n.iter : 10
random.seed : 0
group.singletons : TRUE
verbose : FALSE
...

```

Re: update

I think it depends on how you built the object. `FindNeighbors.SCT.pca` should be result of running `sctransform normalization`. If you use the *standard way*, you may have other names (e.g.`FindNeighbors.pca`).
> 1 votes

---
Tags: seurat, single-cell, clustering
---


# Question
Title: pysam or piping samtools view to a python script
Is it faster to use the PySam package to run a python script on a bam

```
for read in samfile.fetch('chr1', 100, 120):
 print read

```

compared to using a pipe and reading line by line from stdin in a python script.

`samtools view example.bam | python script.py.`

```
for line in tqdm(sys.stdin, total=number_of_lines)

```

# Answer
> 5 votes
It's a bit hard to say with certainty, though I would suspect that offloading the BAM decompression by using a pipe will be very slightly faster. Note that decompressing and parsing the BAM file will not be the bottleneck in your processing, rather the python script itself will be.

---
Tags: python, sam, samtools, pysam
---


# Question
Title: Bioinformatics flowchart for Materials and Methods
I am going to start work on the Materials and Methods for a paper I contributed to a lot in the bioinformatics processing part. I have seen some papers nowadays tend to first present the processing steps that they followed as a flow chart, to show the reader what happened generally before going into details in the written M&M section.

Are there any specific guidelines for creating this flowchart? Do you have any tips or advice from personal experience?

---

(I am aware that this may be off-topic. I have a meta question about it here)

# Answer
If you're using a workflow management system for pipelining like snakemake or ruffus, you can actually build a flowchart based on the rules/functions in your code:

Snakemake example with the `--dag` flag:

Ruffus example using `pipeline_printout_graph`:

You'll probably still need edit the image aftwards in software like inkscape or illustrator though.
> 4 votes

---
Tags: literature
---


# Question
Title: Why is it necessary to add hydrogen and delete water before protein-ligand docking?
What is the reason for adding hydrogen and removing unnecessary water molecules from the protein structure before protein-ligand docking? FYI, the tools I used for docking is GOLD.

# Answer
> 4 votes
My answer is mainly based on molecular dynamics,but in this context it shouldn't make much of a difference.

The main reason has to do with the `force-field`-`water` model used during the simulation. Examples of water models are `SPC`, `TIP3P` or `TIP4P`; same molecule (water) but different topology.

A `pdb` usually comes from different sources and it can contain molecules that are not standard and therefore not recognised by the model of your simulation. Usually, once the structure is loaded, your system is `re-built` so that it conforms to the model topology. There are exception though, for example if your water is part of a functional active site. In that case you can keep it, but its structures (topology of the water molecule) has to be recognised by the model.

Hydrogens are usually added to "saturate" and stabilise the bonds.

# Answer
> 4 votes
## Hydrogens

Most crystallographic structures are X-ray and the hydrogens are absent. In the case of neutron diffraction models, there may be many missing high b-factor ones (example in PDB:5MOP). As a consequence models are reprotonated. Also missing carbons and missing loop residues are added for some in silico applications too (not docking).

However, it is always good to manually tweak the protonation. In the case of enzymes with protanatable side chains in the active site, some reading is required to get the charges right: by default aspartates will be unprotonated, but they may be protonated (glutamic acid) in the active site and it is important to get it right (*cf.* this post about it). Some times its the ligand parameterisation may give the physiologically incorrect protonation (*cf.* this post).

## Waters

Crystallographic waters can be important catalytically and structurally, but can be pushed out and actually it is a good thing: you have a hit that take its place. Polar bonds are stronger than hydrophobic ones and most types of π bonding.

Waters in most docking software is `implicit` and is modelled as a parameter in the force-fields —think of it as a ghostly watery haze. So a left behind water, not only will get in the way of the docked small molecule, it may for some programs cause issues as they will try to dock it. A notable exception is WaterDock mod for AutoDock, but the jury is out on whether it is better. Water is polarisable in reality in a very complex way (protonated (OH3+) or deprotonated (OH-) with very different charges) and most water models do not cut it and instead do not account for this. The way implicit waters are computed differs in some applications, which is also something to look into.

# Answer
> 0 votes
The molecular docking programs usually tread the target structure as rigid. It does not know which atom belongs to the receptor and which is bulk water. So, if you keep the water in the structure the program cannot properly place the ligand there.

---
Tags: protein-structure, docking, structural-biology
---


# Question
Title: bash loop for strelka2 processing for multiple files from one directory
In the below `bash loop for strelka2`, 3 bam files are being stored in a variable. The two questions are:

Is this proper way to loop through multi bam and wait for workflow to be created for `strelka2` processing, and can you create a workflow in a directory without the bam in it. Thank you :).

**contents of dir**

```
s1.bam
s2.bam
s3.bam

```

**bash**

```
dir=/path/to/bam
data=/path/to/processed/data

cd $dir
 for bam in ${dir}/*.bam; do sample_basename=$(basename "${bam}"); done # loop through dir and get bam

${path_to_strelka}/bin/configureStrelkaGermlineWorkflow.py \
  --bam ${dir}/${sample_basename} \
  --referenceFasta <fasta> \
  --runDir ${data}
  wait  # wait until workflow in data is created before executing below  in data
${data}/runWorkflow.py -m local -j 20

```

# Answer
> 1 votes
The loop you're using does not work, as it overwrites the value of `sample_basename` each time. Plus, given that you're passing in the full path to the bam (using `$dir/$sample_basename`) anyway, I don't see the reason to use the `sample_basename` variable at all.

Just this:

```
dir=/path/to/bam
data=/path/to/processed/data
ref_fa=/path/to/ref/fasta

cd $dir

for bam_file in *.bam
do
  ${path_to_strelka}/bin/configureStrelkaGermlineWorkflow.py \
    --bam ${bam_file} \
    --referenceFasta ${ref_fa} \
    --runDir ${data}
  ${data}/runWorkflow.py -m local -j 20 #I'm assuming this goes inside the loop. If not, move it outside the loop. You don't need a `wait` in a shell script that's working line-by-line already.
done

```

---
Tags: bash
---


title: Density plot, scale it to 0-1

text:
# Question
Title: Density plot, scale it to 0-1
Doubt regarding density plot what is the scale being plotted in the Y axis.Is it possible to scale it to 1?.

This is my code

```
w <- read.csv("Normal_Myeloid_Dev_stages/Myeloid_non_coding_non_CDS_CORRVALUE.txt",header = TRUE)
head(w)
#names(w)[1] = "Sample"
head(w)

#w.plot <- melt(w)

df <- w

df_melt=melt(df,id.vars="Sample")
#head(df_melt)
tail(df_melt)
df_melt$Group <- gsub('[0-9]', '',df_melt$variable)
head(df_melt)

p1 <- ggplot(aes(x=value, colour=Group), data=df_melt)
p1 + geom_density()

ggplot(df_melt, aes(x = value,y=..density.., fill = Group)) + geom_density(alpha = 0.5)+
  theme_bw(base_size=30)+
  theme(axis.text.x=element_text(angle = 45, size=45, face="bold", hjust = 1),
        axis.text.y=element_text(angle=0, size=50, face="bold", vjust=0.5),
        plot.title = element_text(size=40, face="bold"),
        legend.title=element_blank(),
        legend.key.size=unit(1, "cm"),      #Sets overall area/size of the legend
        legend.text=element_text(size=30))

dev.off()

```

Second question regarding finding the overlapping area between various condition let say I have four condition and i want to find out what is the overlapping or intersecting area between them .

I tried on of the stack solution firstly is it the right way to do it for finding overlap? if yes can i try finding all the condition overlap such as I have HSC,CMP,GMP and Mono.

```
FDensity = approxfun(density(df_melt$value[df_melt$Group=="CMP"], from=.80, to=1))
MDensity = approxfun(density(df_melt$value[df_melt$Group=="Mono"], from=.80, to=1))
plot(FDensity, xlim=c(.80,1), ylab="Density")
curve(MDensity, add=TRUE)

FminusM = function(x) { FDensity(x) - MDensity(x) }
Intersect = uniroot(FminusM, c(.8, 1))$root
points(Intersect, FDensity(Intersect), pch=20, col="red")

integrate(MDensity, .80,Intersect)$value +
  integrate(FDensity, Intersect, 1)$value

```

Any suggestion or help would be really appreciated

My dataframe small subset

```
dput(head(df))
structure(list(Sample = structure(c(9L, 10L, 11L, 12L, 1L, 2L
), .Label = c("CMP1", "CMP2", "CMP3", "CMP4", "GMP1", "GMP2",
"GMP3", "GMP4", "HSC1", "HSC2", "HSC3", "HSC4", "Mono1", "Mono2",
"Mono3", "Mono4"), class = "factor"), HSC1 = c(1, 0.901758000194052,
0.880971458703505, 0.900098712568466, 0.873507054094022, 0.906315987662777
), HSC2 = c(0.901758000194052, 1, 0.945122894369186, 0.955453795442752,
0.844777154811663, 0.960260413763721), HSC3 = c(0.880971458703505,
0.945122894369186, 1, 0.931891528389177, 0.829521543809914, 0.93636611655036
), HSC4 = c(0.900098712568466, 0.955453795442752, 0.931891528389177,
1, 0.852155096935692, 0.947192034704188), CMP1 = c(0.873507054094022,
0.844777154811663, 0.829521543809914, 0.852155096935692, 1, 0.88406635204624
), CMP2 = c(0.906315987662777, 0.960260413763721, 0.93636611655036,
0.947192034704188, 0.88406635204624, 1), CMP3 = c(0.883085447599108,
0.928929790213059, 0.937278860333014, 0.923139729112196, 0.870595169861622,
0.958330253497026), CMP4 = c(0.903500691840647, 0.931015449611016,
0.915510499719426, 0.949501343959892, 0.89151869084162, 0.960206335031192
), GMP1 = c(0.859320159223793, 0.833454007185579, 0.821479398953591,
0.837605565046685, 0.856892209387879, 0.8612685942226), GMP2 = c(0.894680165774456,
0.944129985346879, 0.92461599550388, 0.937108867595306, 0.875930650641323,
0.966512850930959), GMP3 = c(0.881416749993221, 0.925001567707159,
0.923015420025178, 0.920758302664238, 0.863998116268766, 0.95069438729762
), GMP4 = c(0.887230702792805, 0.911614382740949, 0.898704725226052,
0.938310378760322, 0.875956580458754, 0.93865124308053), Mono1 = c(0.80834274546097,
0.803706049167148, 0.785535924609134, 0.794438848049474, 0.794402126200663,
0.810382054101301), Mono2 = c(0.829928035079747, 0.868021129791583,
0.848929312833536, 0.852953184941131, 0.818074585324121, 0.880490870699251
), Mono3 = c(0.834028389107401, 0.859818911060746, 0.852694253707916,
0.85039434449057, 0.82251614206182, 0.872740655210769), Mono4 = c(0.832524159512094,
0.858444209302146, 0.84397886802844, 0.862850406374668, 0.818436260136434,
0.870293804485003)), row.names = c(NA, 6L), class = "data.frame")

```

# Answer
> 6 votes
I'm not fully to have understood your second questions, but first, if you want to scale the density to 1, you can use `y= ..scaled..` in your `aes`:

```
library(tidyverse)
df %>% pivot_longer(., everything(), names_to = "Variable", values_to = "Value") %>%
  ggplot(., aes(x = Value, y = ..scaled.., fill = Variable)) + geom_density(alpha = 0.4)

```

Now regarding the second question, if you are looking for the overlap between each of your density, you can use the package `overlapping` as specified in the post you are referring: https://stackoverflow.com/questions/41914257/calculate-area-of-overlapping-density-plot-by-ggplot-using-r?noredirect=1&lq=1

```
library(overlapping)
library(lattice)
x <- list(D1 = D1, D2 = D2, D3 = D3, D4 = D4)
overlap(x, plot = TRUE)

```

And if you want to extract the value of overlap, you can get it:

```
OV <- overlap(x)$OV

> OV
    D1-D2     D1-D3     D1-D4     D2-D3     D2-D4     D3-D4
0.5971760 0.3171644 0.7566528 0.4926057 0.6683083 0.3756224

```

Is it what you are looking for ? Or did I misunderstood your question ?

**Dummy Dataset**

```
# Dummy example:
D1 = rnorm(100, mean =1, sd = 2)
D2 = rnorm(100, mean = 2, sd = 2)
D3 = rnorm(100, mean = 3, sd = 1)
D4 = rnorm(100,mean = 1.5, sd = 1.5)

df = data.frame(D1,D2,D3,D4)

```

---

**EDIT: Using reproducible example provided in the question** For plotting `geom_density` of your data:

```
library(tidyverse)
df %>% pivot_longer(., -Sample, names_to = "Variable", values_to = "Value") %>%
  mutate(New_Var = gsub("\\d","",Variable)) %>%
  ggplot(., aes(x = Value, y = ..scaled.., fill = New_Var))+
  geom_density(alpha = 0.4)

```

And to find the overlap:

```
d<- df %>% pivot_longer(., -Sample, names_to = "Variable", values_to = "Value") %>% mutate(New_Var = gsub("\\d","",Variable))
a <- split(d$Value, unique(d$New_Var))

g <- overlap(a, plot = TRUE)$OV

```

I don't think there is an option in the `overlap` function to scaled each overlap plot to 1 but these plot are using `ggplot` and so, I imagine you can extract their points and scaled them by plotting them separately in `ggplot2`.

Hope it answers your question

---
Tags: r, ggplot2
---


# Question
Title: Finding simple sequence from reads with significant overlap
I wrote a "script" to pull out reads from a huge fastq file in an iterative manner, by finding homology to the previous sequence.

It should be relatively easy to overlap them and assemble the gene, however traditional assemblers like SPAdes fail due to *obviously* low coverage.

I am hesitant to just feed it more faked reads, when I know I have all the data necessary to just overlap them into a single sequence.

The reads are set up like this:

```
GATCAAACATC
     AACATCAGTTAG
         TCAGTTAG
              TAGAGGATAGC

```

What kind of assembler is out there where I can feed it this kind of terribly-formatted data?

# Answer
Miniasm will probably work for this; it's designed to be an ultra-fast, greedy assembler, which has been used for people doing real-time assemblies from MinION data. Because of its greedy nature, it can misassemble regions where the overlapping fragments are ambiguous, so it's worth using some other method in addition to miniasm to confirm that the assembly is correct.

By default it looks for overlap regions that have a coverage of 3 reads when joining (option `-c`), but I assume that could be reduced down to 2 reads, as in your example (which has a minimum overlapping coverage of 2).
> 1 votes

---
Tags: sequence-alignment, assembly
---


# Question
Title: What score value is considered acceptable when using Vina Autodock tool for molecular docking simulations?
What score value is considered acceptable when using Vina Autodock tool for molecular docking simulations?

# Answer
> 1 votes
The answer is "Depends". But the best thing is always to look at what you ended up with and do some triaging before buying/synthesising anything. By triaging I mean filtering by its properties (*.e.g.* Lipinski, reactivity, synthetic feasibility *etc.* —discussed in this Bioinfo SE question).

## Score

Regardless of docking software used, a forcefield free-energy score is a result of the scoring function used and of the properties of the molecule —including VENA which uses a hybrid score function.

This score is *not*, but generally is set up to try to match, the actual phys-chem one. Gibbs free energy is measured in kcal/mol (or less commonly in kJ/mol). A hydrogen bond is about 2 kcal/mol, a π-π bond 1 kcal/mol and a salt bridge 3. The average energy of a collision of RT water molecules is 0.6 kcal/mol.

Therefore, if you dock multiple compounds (virtual library screen, VLS) you sort by free energy to see the strongest binders and consider the top ones. However, a catch is that this will be biased by MW.

Also, it is useful to score the hits with a different program (say Gold or ICM-dock).

## Ligand Efficiency

A variant metric that is employed when your hits vary in size and when they are over a given size (*e.g.* not fragments under 200 Da) is Ligand Efficiency (LE), which is the ratio of your score over the number of heavy atoms (ref).

## Single compound

If you have a single compound and want to know if it binds acceptably, then this is a length-of-a-piece-of-string question as every compound has an affinity constant to a given protein, except most may be beyond their solubility. Initial hits will have µM affinity, while a lead will have nM affinity. However, a KD is not a ∆G and depends on your protein.

---
Tags: docking
---


# Question
Title: Merging files gives memory allocation error
I have two files like

```
> head(a)
  chr position
1   1   136962
2   1   562020
3   1   672948
> class(a)
[1] "data.frame"
>

```

And

```
> dim(b)
[1] 1855235      10
> head(b[,1:3])
            SNP chr position
1:  rs62635286          1    13116
2:  rs75454623          1    14930
3:    rs806731          1    30923
4: rs200943160          1    49298
5: rs116400033          1    51479
6: rs141149254          1    54490
> class(b)
[1] "data.table" "data.frame"
>

```

I want to find corresponding SNP in file b matched to position in file a so I did like this that gives memory allocation error.

```
> unique(merge(a,b,x.by=a$position,y.by=b$position))
Error: cannot allocate vector of size 137.2 Gb

```

I also tried rstudio cloud but uploading the bigger file never completes

Can you help me with task please?

The expected output would be something like

```
rs62635286  1   13116
rs376723915 1   91515
rs147061536 1   92858
rs62642131  1   135982
rs371474651 1   158006
rs201293782 1   665401
rs138476838 1   668374
rs1401137   1   2063094
rs3128293   1   2065339

```

> EDITED

```
> setDT(a)
> setkey(a, chr, position)
> setkey(b, chr, position)
Error in setkeyv(x, cols, verbose = verbose, physical = physical) :
  some columns are not in the data.table: chr,position
> setkey(b, Chromosome, Position)
> a[b, nomatch=0]
Error in bmerge(i, x, leftcols, rightcols, roll, rollends, nomatch, mult,  :
  Incompatible join types: x.chr (factor) and i.Chromosome (integer). Factor columns must join to factor or character columns.
>

```

# Answer
> 2 votes
You could try something like this :

First you should work with the same class of object : `data.table` for fast merging operation

```
# Turning a to a data.table
setDT(a)
# b is already a data.table

```

Then you might have to set the keys of both table :

```
# I added chr in case you need to match the same position on a different chromosome.
setkey(a, chr, position)
setkey(b, chr, position)

```

Finally perform the joining operation :

```
a[b, nomatch=0]

```

Because you didn't provide the expected output, it might need some tuning !

---
Tags: r
---


# Question
Title: Generate SNP/indel annotation in Arabidopsis
I have a genome-wide list of germline SNPs and short indels for Arabidopsis thaliana, which I generated with Varscan. Regardless of the tool used to generate them, I would like to annotate them, i.e. knowing which ones can cause an aminoacid change or an early stop, using the default Arabidopsis thaliana Columbia 0 cultivar, for which I have both the sequence (FASTA from TAIR10) and the updated annotation (GFF from Araport).

```
Chrom   Position    Ref Var
Chr1    626503  G   T
Chr1    926694  C   T
Chr1    5280350 C   A
Chr1    5699993 C   A
Chr1    7004559 G   A
Chr1    8325810 C   T
Chr1    9371723 T   G

```

What I want to do is similar to what Annovar does, but unfortunately Annovar does not support Arabidopsis. I was thinking of an already existing R pipeline that takes in a genome, an annotation, a SNP/indel list and *boom*, annotation. But I couldn't find any, except maybe snpEffect. Any tips?

# Answer
> 2 votes
The Ensembl Variant Effect Predictor will annotate your A. thaliana variants: http://plants.ensembl.org/Oryza\_sativa/Tools/VEP

The web based tool will use the reference genome assembly and gene annotation used in the Ensembl browser: http://plants.ensembl.org/Arabidopsis\_thaliana/Info/Annotation/#assembly

However, if you use the command-line version of VEP, you can provide your own reference genome assembly and gene annotation files to be used as a reference.

You can find all of the documentation here: http://www.ensembl.org/info/index.html

Best wishes

Ben Ensembl Helpdesk

# Answer
> 1 votes
I recently used SIFT4G to predict the deleteriousness for a SNP dataset in Arabidopsis.

I noticed that besides this, the tool also provides annotation for the SNPs (like gene region, amino acid substitution, start/stop etc.)

They have an Arabidopsis database readily available so you don't have to generate your own. Your input should be enough to run the tool, although you might have to do some reformatting (the tool requires VCF format files).

One thing though, not all SNPs are considered by the tool. Intronic SNPs and SNPs outside gene regions are ignored afaik.

---
Tags: sequence-annotation, snp, indel
---


# Question
Title: how to change the PC use in the dimplot and feature plot
I would like to know how to change the PC use in the dimplot and featureplot by using Seurat. for we can get the x-axis and the y-axis like PC-1 and PC-2, if I want to use PC-4 and PC-5. How to modify the code?

# Answer
You can choose which PC dims to plot by specifying the `dims` argument within the functions. To use PC-4 and PC-5,

```
DimPlot(object = pbmc_small, dims = c(4, 5), reduction = "pca")

```

and

```
FeaturePlot(object = pbmc_small, features = "Xist", dims = c(4, 5), reduction = "pca")

```

---

Seurat documentation is available here.
> 1 votes

---
Tags: seurat, single-cell
---


# Question
Title: nanopore - where to retrieve information from the basecaller used
Where can I find information about the basecaller used in a specific nanopore run (we have a Gridion machine)?

I know already that we have `Guppy v. 3.0.3`, which I found out by clicking on the `about` tab on the MinKNOW GUI).

I was wondering, though, if additional information about the basecaller is outputted in any of the files generated when the run is finished? I checked already the `final_summary.txt` and it does not contain such information. I also checked the head of the `fastq` files and no further details.

My ultimate goal for the moment is to run medaka on a draft genome assembly. But for that I'd need to specify the `medaka model`. From `medaka_consensus -h` I see that the options available and defaults are:

```
-m  medaka model, (default: r941_min_high_g330).
        Available: r941_min_fast_g303, r941_min_high_g303, r941_min_high_g330, r941_prom_fast_g303, r941_prom_high_g303, r941_prom_high_g330, r10_min_high_g303, r10_min_high_g340, r941_prom_diploid_snp, r941_min_high_g340_rle.
        Alternatively a .hdf file from 'medaka train'.

```

Therefore I already know that my option should contain `min` and `g303` (format of the models available described here). But know I'd need to find out the `{pore}` and `{caller variant}`

# Answer
I basecall separately with guppy. In the output folder specified by `--save_path` or `-s` there are a whole bunch of .log files.

```
$ ls -l *.log | head
-rw-r--r-- 1 tom tom 5242714 Dec  3 11:04 guppy_basecaller_log-2019-12-02_22-02-36.log
-rw-r--r-- 1 tom tom 5242718 Dec  3 11:06 guppy_basecaller_log-2019-12-02_22-04-38.log
-rw-r--r-- 1 tom tom 5242730 Dec  3 11:08 guppy_basecaller_log-2019-12-02_22-06-41.log
-rw-r--r-- 1 tom tom 5242879 Dec  3 11:11 guppy_basecaller_log-2019-12-02_22-08-50.log
-rw-r--r-- 1 tom tom 5242811 Dec  3 11:13 guppy_basecaller_log-2019-12-02_22-11-06.log
-rw-r--r-- 1 tom tom 5242817 Dec  3 11:15 guppy_basecaller_log-2019-12-02_22-13-28.log
-rw-r--r-- 1 tom tom 5242880 Dec  3 11:18 guppy_basecaller_log-2019-12-02_22-15-47.log
-rw-r--r-- 1 tom tom 5242816 Dec  3 11:20 guppy_basecaller_log-2019-12-02_22-18-09.log
-rw-r--r-- 1 tom tom 5242777 Dec  3 11:22 guppy_basecaller_log-2019-12-02_22-20-30.log
-rw-r--r-- 1 tom tom 5242742 Dec  3 11:25 guppy_basecaller_log-2019-12-02_22-22-47.log

```

The information you need for `medaka` is at the top of the **first** log file (alphabetically, and by creation date).

```
$ head -n 13 guppy_basecaller_log-2019-12-02_22-02-36.log
2019-12-02 22:02:36.255793 [guppy/message] ONT Guppy basecalling software version 3.4.1+ad4f8b9
config file:        /opt/ont/guppy/data/dna_r9.4.1_450bps_hac.cfg
model file:         /opt/ont/guppy/data/template_r9.4.1_450bps_hac.jsn
input path:         raw/minion
save path:          test_341/minion
chunk size:         1000
chunks per runner:  512
minimum qscore:     7
records per file:   4000
num basecallers:    1
gpu device:         cuda:0
kernel path:
runners per device: 4

```

If you ran live basecalling and you don't have those log files you could also try a recursive grep from the output folder, e.g.

```
grep -nr "^model file" .

```
> 2 votes

---
Tags: nanopore, assembly, long-reads, base-calling
---


# Question
Title: Given a .gb file and a locus - how to get relevant annotations in Python?
Given a `.gb` file and a specific locus in the genome - how can I parse the relevant annotations in Python (i.e., annotations that include that locus)?

I import the features using:

```
SeqIO.read('my_gb_file.gb', 'gb').features

```

I then wish scan/filter sequence meta data.

Is there a function or method in `Biopython` that does this?

# Answer
> 1 votes
Biopython is the main package for this. It's only a few lines, so it is not reinventing the wheel.

So you want to iterate across the features table (a list in Biopython) of a record and find the case where the qualifier\['locus'\]\[0\] matches your query.

The things to watch out for are:

* filter by `type` (CDS?)
* entries where there is not key among the qualifiers, so add error catching (try: ... except Exception: pass
* the values of the qualifiers of a feature are a list. so add a \[0\]
* multithreading.Pool might help if speed is a worry (which it really shouldn't). Asyncio and thread are not the way (both single core).

A nice handy way to do operations on the `features` is to convert it to a pandas dataframe.

```
from Bio.SeqFeature import SeqFeature
from Bio import SeqIO
from Bio.Seq import Seq
import pandas as pd

        def convert(feat: SeqFeature, flatten:bool=True):
    if flatten:
        clean = lambda v: v[0] if isinstance(v, list) else v
        qualifiers = {k: clean(v) for k, v in feat.qualifiers.items()}
    else:
        qualifiers = feat.qualifiers
    return {'id': feat.id,
            'type': feat.type,
            'location_strand': feat.location.strand,
            'location_start': feat.location.start,
            'location_end': feat.location.end,
            **qualifiers}

transcript = SeqIO.read('.gb', 'genbank')
feats = pd.DataFrame([convert(feat) for feat in transcript.features])

```

This means that in a jupyter notebook one can display the table to look at it (`feats` as last line in a cell or `IPython.display.display(feats)`), export it to a csv (`feats.to_cvs('.csv')`) or subset it:

```
subfeats = feats.loc[(feats['type'] == 'exon') & (feats['location_start']  > 100)& (feats['location_end']  < 500)]

```

---
Tags: sequence-annotation, python
---


# Question
Title: Error in apply(dat[, batch == batch_level], 1, function(x) { : dim(X) must have a positive length
I am new here so forgive me if I don't provide all the details in the first go. There is a very similar question (how to solve “dim(X) must have a positive length” at running ComBat function in R) but the solution did not work for me.

I am getting `Error in apply(dat[, batch == batch_level], 1, function(x) { : dim(X) must have a positive length` when I use the `sva` library's `comBat()` function.

It works when I use the following:

```
> dim(log_proteins_norm_sorted)
[1] 5415   66

> log_proteins_norm_sorted[1:5,1:5]
      1027     1075     1116     1119    11211
1 22.50513 20.64312 21.56182 22.01455 23.17497
2 21.67946 20.90994 21.61928 21.92025 21.98503
3 21.06812 20.52647 20.72076 20.79732 21.38957
4 18.95485 19.70027 18.39504 20.20218 19.79887
5 22.75966 22.32306 22.32936 22.39271 22.84204

> dim(samples_pheno)
[1] 66 24

> samples_pheno[1:3,1:3]
     Trial Number  Trial Sample
1027      222       A   1027
1075      223       A   1075
1116      224       A   1116

> proteins_combat = ComBat(dat= as.matrix(log_proteins_norm_sorted), batch=samples_pheno$Batch ,  par.prior = T)
Found10batches
Adjusting for0covariate(s) or covariate level(s)
Standardizing Data across genes
Fitting L/S model and finding priors
Finding parametric adjustments
Adjusting the Data

```

So as you can see, this works. However, when I subset some of the samples and reduce the sample number from 66 to 55 and run comBat, I get the following:

```
> dim(log_proteins_MRD_sorted)
[1] 5415   55
> log_proteins_MRD_sorted[1:5,1:5]
      1027     1075     1116     1119     1217
1 22.55492 20.69292 21.61162 22.06435 21.86900
2 21.72926 20.95974 21.66907 21.97005 21.35203
3 21.11791 20.57627 20.77055 20.84712 21.01896
4 19.00465 19.75006 18.44483 20.25198 18.95107
5 22.80946 22.37285 22.37915 22.44250 22.50524

> dim(samples_pheno3)
[1] 55 24
> samples_pheno3[1:3,1:3]
     Trial Number  Trial Sample
1027      222        A   1027
1075      223        A   1075
1116      224        A   1116

> proteins_combat_MRD = ComBat(dat= as.matrix(log_proteins_MRD_sorted), batch=samples_pheno3$Batch ,  par.prior = T)
Error in apply(dat[, batch == batch_level], 1, function(x) { :
  dim(X) must have a positive length

```

In case you wanted to know, I am working with the same original file for `log_proteins_MRD_sorted` and `log_proteins_norm_sorted`. I just removed the 11 samples from the 66 samples and then normalised the data again (hence slight difference in values between the two objects.

For `samples_pheno3`, I run this code:

```
> #Subset the phenotype data to those samples with only MRD data
> samples_pheno2 <- samples_pheno[-c(5,6,11,16,21, 26, 29,30,32, 45,56),]
> MRD_NA1 <- rownames(samples_pheno2)[!is.na(samples_pheno2$MRD)]
> samples_pheno3 <- samples_pheno2[MRD_NA1,]

```

Also this is out my data looks like for the batches:

```
> samples_pheno3$Batch
 [1] 7  5  6  9  7  7  7  7  8  8  8  8  9  9  9  9  10 10 10 10 11 11 11 2  7  2  2  2  9  3  3  3  3  3  3  3  7
[38] 10 3  4  10 4  4  8  4  4  4  5  5  5  8  11 5  5  5
Levels: 2 3 4 5 6 7 8 9 10 11

```

So just by subsetting the data to 55 samples, it doesn't seem to work for some reason. Is it the way I am subsetting the data which is confusing the code?

I would be grateful for your advice.

# Answer
> 1 votes
With the help of my very capable local bioinformatician, I have found a solution. Essentially, the problem was that when I reduced the sample size from 66 to 55, there was a batch number which only had 1 sample in it. When I removed that sample from the analysis, and ensured that each batch had at least 2 samples in it, the `ComBat()` then worked nicely!

Hope this helps future googlers.

---
Tags: combat, sva
---


# Question
Title: Parallelizing WGCNA k-means clustering and merging smaller clusters
WGCNA carefully takes advantage of server parallelization in several functions, and I can switch this on with:

```
enableWGCNAThreads(nThreads=28)

```

That speeds up some functions significantly, but the two bottleneck functions in the pipeline (in my experience) are part of the "Projective k-means" operation, primarily `"..k-means clustering.."` and `"..merging smaller clusters.."`.

Unfortunately, parallelization does not work for those two functions.

Am I missing something? How can I parallelize these functions? Otherwise I'm wasting HPC resources and cluster time.

For reference, I'm working on a 28-processor HPC node with 96GB RAM, and I'm trying to run analysis on a 15,000-square correlation matrix.

# Answer
The bad news is that, indeed, projectiveKMeans is not parallelized and I am not sure how much of it is (easily) parallelizable. The good news is that with 15k features (genes) and 96GB RAM you don't need to run preclustering at all. Just analyze the whole data set in one block.

If you for some reason insist on splitting the data into blocks, make sure the BLAS R uses on the cluster is a fast BLAS, ideally one that can parallelize matrix operations. A lot of time in merging of small clusters is spent calculating singular value decompositions, which becomes much faster using a good BLAS. HPC clusters will often have Intel KML available, or the admins could install OpenBLAS.
> 2 votes

---
Tags: r, wgcna
---


title: What is a read count?

text:
# Question
Title: What is a read count?
A very basic question: What exactly is a read count in shotgun sequencing and how is it calculated? I'm particularly interested in understanding the definition of the read count in the context of shotgun metagenomic sequencing though I assume the same concept applies to other types of shotgun sequencing as well. Surprisingly, a Google search did not results in any answers!

In the context of metagenomics, the taxonomic profiling of metagenomic data provides a list of microbial taxa and their "counts". These counts are often normalized to account for genome size or are normalized to convert counts to proportions. Understanding the concept of a read count may help me better understand all these. Thanks!

# Answer
To add to the answer from Bioathlete, "reads" in the traditional genomic sense are the string of nucleotide data inferred from the sequencer from a given fragment of genetic material (could be Illumina, Nanopore, PacBio, 10X, etc.).

So read count first applies for all of the "reads" defined above generated by the sequencer. From a NovaSeq, the total number of bases from a single flowcell is up to 3 trillion bases at 150bp each read for the highest throughput flowcell (so a read count of 3,000,000,000,000/150 = 20B total read count).

The read count for a single Nanopore flowcell is roughly (30 billion bases/avg read length), but the read lengths can be SUPER long (100kb compared to 150bp is substantial). Long read technologies are generally used for different analyses than short reads.

So all of the above reads would be defined as "raw reads", meaning they are reads not yet analysed in any way to be used for data analysis. So using illumina as an example, for a single Whole Genome sequencing sample we might have 800,000,000 raw reads. We then use an aligner like bwa mem to map the reads to a single reference genome (or in your case you would use a different aligner to map to multiple genomes). From here we still have the concept of raw reads, but we also now have the total number of mapped reads. A typical good illumina library for a single reference genome will have \>90% mapping rate for a human sample.

We can then say we have 800,000,000 raw reads, and 720,000,000 mapped reads with a mapping to a reference genome of 90%. So we can compare "raw read count" to "mapped read count" to get an idea of the quality of our library.

In the case of metagenomics, you would consider your "raw read count" to still be the total number of reads you got off the sequencer, but the "mapped read count" would be the sum of the reads mapped to each individual reference genome to get a total quantity of "mapped reads" which could give you an indication of library quality. I'm not sure what a good metagenomics library "mapped read count" is, but I'm sure you can find that information.

Lastly, we want to get an idea of the depth of coverage for each mapped read to the genome. This will tell us how well our mapped reads actually covered the reference genome or genomes we are interested in.

I hope this was helpful!
> 5 votes

# Answer
Typically read count is the total number of reads going into the analysis. It could be based off single or multiple sequencing libraries.

Also it can be used to describe the number of reads that align to a region of the reference. Depth or coverage are also terms used in this case.

You will need to be more specific in your use case for an accurate definition.
> 2 votes

---
Tags: ngs
---


title: Comparing factor levels in deseq2

text:
# Question
Title: Comparing factor levels in deseq2
I have sample information (sampleinfo) in this format:

```
SampleID    condition   type
MP_1    A_1 Paired-End
MP_2    A_1 Paired-End
MP_3    A_1 Paired-End
MN_1    A_2 Paired-End
MN_2    A_2 Paired-End
MN_3    A_2 Paired-End
WP_1    B_1 Paired-End
WP_2    B_1 Paired-End
WP_3    B_1 Paired-End
WN_1    B_2 Paired-End
WN_2    B_2 Paired-End
WN_3    B_2 Paired-End

```

I want to compare A\_1 with A\_2, B\_1 with B\_2 ans also want to compare (A\_1+A\_2) with (B\_1+B\_2).

This is how I used DESeqDataSetFromMatrix

```
ddsMat <- DESeqDataSetFromMatrix(countData = seqdata[,rownames(sampleinfo)], colData =  sampleinfo, design = ~ condition)
ddsMat <- DESeq(ddsMat)

res.05_FC0<-results(ddsMat,contrast = c("condition","A_1","A_2"),alpha=0.05,lfcThreshold = 0)
res.05_FC0<-results(ddsMat,contrast = c("condition","B_1","B_2"),alpha=0.05,lfcThreshold = 0)

```

How can I compare (A\_1+A\_2) with (B\_1+B\_2) using the contrast argument in results function?

Of course I can make a new sampleinfo table, but I want to use the same sampleinfo as above. I tried to read about numeric contrast, but I am unable to use it correctly.

# Answer
> 1 votes
You can use `contrast = c(0.5, 0.5, -0.5, -0.5)`, which means "the average of the first two coefficients minus the average of the second two.

---
Tags: deseq2
---


title: Use of Electronic Phenotype in EHR

text:
# Question
Title: Use of Electronic Phenotype in EHR
May I know what's the use of `Electronic Phenotyping using EHR data`?

I did refer this link but have few questions

I understand that `Phenotype` is a set of criteria that you apply on your EHR data to select patients of interest from EHR database.

However what I am trying to understand is the difference between ordinary SQL query and electronic Phenotype

If both are same,

1) what I would like to know is the use of `electronic phenotyping`? How it is used?

2) Why do we have `pheKB` repository and other repositories? What's the use in validating them in multiple sites?

3) Is it like identifying patients real time? I am not able to get the usefulness of this. I understand that we select patients of our interest (cohort), but i can just write an ordinary SQL and send it to my EHR database to identify patients. but what's special about Electronic Phenotype?

4) I also see that they use multiple conditions (Lab tests, Medications, Diagnosis etc) to identify patients. What am I missing here?

Can someone help me by giving an simple example? I am new to healthcare and trying to learn. Would really be helpful

# Answer
> 1 votes
> I understand that Phenotype is a set of criteria that you apply on your EHR data to select patients of interest from EHR database.

No, a phenotype is a way of behavior or other observed characteristic of a person, resulting from a combination on its underlying biology and the environment it's in. Think things like, height, running speed, strength, etc.

In the context in which you're interested, this would be more properly, "the clinical presentation of a disease or disorder".

1. Step one of diagnosing a patient is defining exactly how they present (i.e., exactly what their maladies, aberrant electrolyte levels, genetic mutations, etc. are).
2. An algorithm that only works on the test population is useless. Google "model overfitting".
3. The point of an algorithm is to take the EHR data and progress to a diagnosis. This isn't real time, you'd just run a batch of algorithms as test results come in.
4. You need patients in your test/validation sets to have the required inputs for your algorithm.

FYI, I don't know of any regulars on this site that deal with EHR data. If there's a community of such people on the web they must be elsewhere.

---
Tags: sequence-alignment, database, genome, statistics, machine-learning
---


# Question
Title: How could I match atom orders between a .mol2 and a .pdb?
When I was drawing conformation of Autodock result, I noticed that Autodock generated pdbqt from dlg has some structural problem (e.g. Benzene ring missing), and wants to correct this. However I noticed that the atom order of pdbqt doesn't match with the original mol2's, and caused some problems when I tried to fix.

So how could I reorder pdbqt according to mol2's atom coordinate?

# Answer
Atom reordering is a common problem in compchemistry.

Rdkit (a python package) can do this, but it is limited by the formats it can read and mol2 files are a bit hit or miss. It works really well with SMILES, SMARTS and mol (sdf) files. But the writing may cause problems with Brookhaven pdb and mol2 files. So the formats you have are both problematic. One is a pdbqt (a PDB file with partial charges (Q, generally Gasteiger calculated) and atom type (T) and the other a mol2. Let's pretend first that they are not.

Also the hydrogens might be made implicit depending on your mol2 file.

```
from rdkit import Chem
ref = Chem.MolFromMol2File('ligand.mol2')
target = Chem.MolFromPDBFile('ligand.pdb')
atoms = target.GetSubstructMatches(ref) #: List[int]
new = Chem.RenumberAtoms(target, atoms)
Chem.MolToMolFile(new, 'neoligand.mol')

```

Extra steps would be required for the protons and if your starting ligand was not from a file, but a SMILES.

```
mol = Chem.AddHs(mol) #protonate explicitly
Chem.GetSSSR(mol) #resonance fixing
AllChem.EmbedMolecule(mol) #initialise for 3d.
AllChem.UFFOptimizeMolecule(mol, maxIters=2000)

```

So this is all nice, but the formats may not work. Luckily openbabel can deal with this problem.

```
obabel -i mol2 input.mol2 -o pdbqt -O output.pdbqt

```

However, one thing to keep an eye out is what will be recalculated. Say, sdf files do not have partial charges or what PDB calls "atom names" and Rdkit calls "atom labels" (About the latter, I actually wrote a script to relabel atom labels in pdbs from mol files), while mol2 and PDBQT have.
> 5 votes

---
Tags: proteins, format-conversion, docking
---


title: Replicate papers (beginner)

text:
# Question
Title: Replicate papers (beginner)
I have a strong mathematical and computer background (also with machine learning) and would like to switch to bioinformatics. I know a fair amount of biology. But I find it hard to find real world practical training in this field. So I would like to try to replicate some bioinformatics papers (with their provided data). Good beginner papers for me would be shorter ones where the biological problem is well stated and the problem is also fitted in to the bigger scope (or is just easily understandable). Do you have any papers to recommend?

# Answer
Bioinformatics is a big and very heterogeneous field with a lot of variation, so it's hard to recommend something that covers everything. That said, I think the best you'll do is to build familiarity with some of the common tools with some canned analyses.

One possible option would be some of the scRNA-seq workflows from the Pachter group, they seem to be fairly conscientious about documentation. Here is one such canned analysis, here is the paper, here is a blog post about it that suggests next steps and discusses at a high level.

scRNA-seq is just what's been trendy for a few years, next year it will be something different. But it at least gets you familiar with some of the file formats and concepts such as alignment, etc.

Hope that helps.
> 1 votes

---
Tags: literature-search, reproducibility, literature
---


title: What is 'k' in sequencing?

text:
# Question
Title: What is 'k' in sequencing?
When a DNA sequence is sequenced, I've only ever dealt with A,T,C,G and N which indicates un-identifiable bases. However, I came across a 'k' recently and I had asked another researcher who gave me an answer for what 'k' represents but I don't quite recall. It can't just be an anomaly in that one sequence file. Any ideas?

# Answer
> 16 votes
See IUPAC codes:

So, as you can see above, `K` means "Either G or T".

# Answer
> 4 votes
It is recommended you learn the degenerate nucleotide code. In sequencing it can signify poor quality sequence data, but in primer design it is useful. R (mutation within purines) and Y (mutation within pyrimidines) are common. K, a purine to pyrimadine or pyrimadine to purine mutation is, in my opinion, rare. I would treat a K mutation with caution and consider the triplet codon around it, i.e. if it is part of a protein gene.

Most phylogeneitcs programs will work with the degenerate nucleotide code, so in theory you can still obtain useful information with it.

---
Tags: sequence-alignment, ngs, phylogeny, phylogenetics
---


# Question
Title: Finding differentially expressed genes between two samples
Hi, I was wondering whether we can find the differentially expressed genes between the `Double-KO` and the `Shox2-KO`. For the `FindMarker()` function, only see the differential expressed gene between different identities. I don't know whether we can find the differential expressed gene between only these two samples.

# Answer
If I understand your question correctly, you can simply use `SetIdent()` to change the "default identity" to your samples and then use `FindMarker()` with the `ident.1 = "Double-KO"` and `ident.2 = "Shox2-KO"`.

Alternatively, you can "collapse" or "summarize" your single cell data into a pseudo-bulk and then use one of the bulk RNA-seq packages.
> 1 votes

---
Tags: seurat, single-cell
---


# Question
Title: How to interpret the GROMACS MD Simulation RMSD Graph?
Hello brilliant researchers! After virtual screening, I performed MD Simulation via GROMACS to check the stability of the protein-top hit complex by following this tutorial http://www.mdtutorials.com/gmx/complex/09\_analysis.html

After that, I facing some difficulties in interpreting the resulted RMSD graph as it is quite different from most of the published RMSD graph.

Based on the graph attached, how to interpret it? Is AChE-dexketoprofen complex (Red) more stable (or conserved) than BuChE-dexketoprofen complex (Black)?

# Answer
The opposite-ish.

Your graph is showing only the ligand, but RMSD is generally used for structures as a whole —so I'll mention about whole protein for future readers.

RMSD is a measure of how much does a structure deviate in position from a reference. Your ligand in the protein in red "wiggles" more than the one in black.

## Structure RSMD

If this were RMSD of the protein (and not the index group for the ligand), it would not really mean that the black protein is more stable. Stability is to do with force-field energetics not position. Some protein have two state positions, say MBP is a flexing dumbbell. This would be seen as a jump in the RSMD graph, similar to that seen for the black trace.

However, the command you used was based on `gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg` from the page you linked. Key here is the index `-n index.ndx`, which limits it to the ligand. Before this happened `gmx trjconv` aligned the protein in the trajectory, so that common mistake was not done.

## Ligand RSMD

In terms of ligand RMSD, lower deviation/less wobble correlates with more stable binding, but energetics (`gmx energy`) will tell you how strongly does it bind.

The black trace has a jump, so your ligand in the black trace (which wobbles less) has an alternative binding mode. This is why you are looking at RSMD. To explore this further you can:

* look at the trajectory
* in PyMOL using the command `split_states` twice, once for the first and then the last state, followed by `align` of the two
* *If it were the whole protein*, looking at the per bond matrix `-bm` option in Gromacs rms and spotting the change.

## 3Å??

I strongly suggest checking the trajectory because both have a deviation of 3Å, which for a ligand is not good.

* Did the energy minimisation converge properly?
* If it is because the proteins themselves wobble substantially, revisit `gmx trjconv` with a focus on the active site minus the ligand.
* Is the ligand binding but leaving one substituent group to wobble?

## Energetics

As with docking, the "score" is free energy difference (in kcal/mol). This is not the same as k\_off and k\_on and their ratio K\_d (Dissociation constant). For those you need to do steered MD or have a way longer simulation (not recommended). The reason being is that the ligand has to make its way into the active site. However, you have a secondary binding mode in black, so have two free energy scores, this means that there is a stronger k\_on.
> 1 votes

---
Tags: proteins
---


title: Adding new columns to this data frame

text:
# Question
Title: Adding new columns to this data frame
I have a dataframe like this

```
> dim(anno_maf)
[1] 32786   129
> head(anno_maf[1:2,1:4])
   Tumor_Sample_Barcode Chromosome Start_Position
1:                    a       chr1       27105930
2:                    a       chr1        3342771
   End_Position
1:     27105930
2:      3342771
>

```

I have a list of genes like

```
gene name
PTEN
PTEN
MAP3K1
MALAT1
MED12
RUNX1
PIK3CA
FBXW7
ARID1A
MED12
TAB1

```

I want to add a column in my dataframe named COSMIC, if a gene in gene list exists in my dataframe ***Hugo\_Symbol*** column, put a **yes** for this gene if not put a ***no***; Finally my data frame has an extra column as ***COSMIC*** in which I have ***yes*** or ***no*** in ***Hugo\_Symbol*** column based on the existence of genes in my gene list

Can you help me please?

Thank you

# Answer
> 2 votes
You can use `ifelse()`:

```
anno_maf$COSMIC <- ifelse(anno_maf$Hugo_Symbol %in% your_gene_list,
                          yes = "yes",
                          no = "no")

```

---
Tags: r, looper
---


title: What does these terms mean?

text:
# Question
Title: What does these terms mean?
I am working on cancer genomics failed to follow some terms even by googling

What is the difference of ***focal genome amplifications*** versus ***convergent amplification*** of given region of genome for instance fibroblast growth factors region?

Thank you for any clarification

# Answer
> 1 votes
I think these definitions still have a degree of subjectivity and I don't think they are mutually exclusive.

Convergent amplification basically means that "similar" cancers/neoplasms have "similar" amplification patterns, i.e. at around the same loci (see Myllykangas S et al.,2006 for an example).

Focal genome amplification is defined as "CNA \[copy number alteration\] of limited size, frequently enriched for cancer driver genes" (Krijgsmanthis et al., 2014).

---
Tags: genomics, wgs, copy-number
---


# Question
Title: Unable to read all alignments in stockholm format : UnicodeDecodeError: 'utf-8'
I have a file in Stockholm format and I am trying to read it using the Biopython package. The code is as follows :

```
import numpy as np
from Bio import AlignIO

alignments = AlignIO.parse("Rfam.seed", "stockholm")
for alignment in alignments:
    print(alignment)

```

It gets read for some alignments and then throws up following error :

```
Traceback (most recent call last):
  File "test.py", line 5, in <module>
    for alignment in alignments:
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Bio/AlignIO/__init__.py", line 394, in parse
    for a in i:
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Bio/AlignIO/StockholmIO.py", line 408, in __next__
    line = handle.readline()
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/codecs.py", line 322, in decode
    (result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xe1 in position 2468: invalid continuation byte

```

How to resolve this error?

# Answer
The last line indicates that this is an encoding problem. Try to see if the function you are using to read the data has an "encoding" parameter that you can set to whatever encoding your file is using. If that doesn't work, you can change the encoding of your file by:

* with `save as` using a text editor of your choice and selecting "UTF8" as the encoding. That would be a quick work-around if you have a small number of files.

or

* You can use the `codecs` library to check and convert file encodings within your code. Here is an example.

Warning: I haven't used `Biopython` in ages but encountered similar problems with Python in the past.
> 2 votes

---
Tags: python
---


# Question
Title: Appropriate value to replace outliers in clinical studies
I am working on clinical study where I am focusing on binary classification problem.

There are few variables which has outliers.

I am trying to know which would be the right value to replace outliers in a clinical study?

Is it `mean`, `median` `min`,`max` or `mode` of a column of interest?

Currently I use `min` and `max` for values below and above 1.5\*IQR threshold

I use `IQR` approach to identify outliers but would like to know how it is done usually in clinical studies?

Can anyone help me with this?

# Answer
In my experience, it is more common to use all *clinical* data as-is for clinical studies. And if data is missing, either omit the sample or omit the variable with missing data.

If your classifier can't handle the wide variation commonly seen in clinical studies then you may want to use a classifier which is less impacted by outliers.
> 3 votes

# Answer
Removing outliers is common practice in statistical modeling and perfectly acceptable. However, with regards 1.5 IQR I am far from certain about this approach. Normally, if you want to be conservative then 3 standard deviations (SD) denote an outlier, which is more stringent than IQR. Some use 2 SD. If the value is lower than 2 SD from the group mean it isn't an outlier, that's because 1.96 SD is the normal distribution.

Normally you would delete the outlier and it takes no further part in the calculation, assigning it to 'na' should work. The subsequent calculation will recognise an 'na' value, I think in python it is np('na') ... Np = numpy.

In this case you appear to be deleting the value and treating it as missing data. Treatment of missing data is complicated I agree 100%. Mean, median are possibilities for missing data I can't see this working for binary data. One approach is regression and in python there is a command specifically to do this (which I forget), R will have the same command and it will assign a value for the missing data. However, you can't do this for outliers.

The other routes to remove outliers are via transformation, log transformation is the most common. Transformation is a big subject.
> 1 votes

---
Tags: public-databases, literature-search, data-management, data-preprocessing, public-dataset
---


# Question
Title: What is the significance of enrichment factor regarding machine learning methods?
I'm fairly new to the field of bioinformatics and ran into a question while reading a paper I found on bioRxiv. The overall setting of the paper is using multi-task deep neural networks for kinase classification. The paper that I'm referring to is *Kinome-Wide Activity Classification of Small Molecules by Deep Learning (Allen et al., 2019)* for anyone who's curious.

I'll get straight into the specific figure that I'm having trouble understanding.

I know that questions like these are a bit discouraged on SE, but I'm having trouble on how to start regarding the interpretation of this figure.

Based on my understanding, the point seems to be that multitaks DNN's have the highest overall EF/EFmax ratio amounts, but what does this ratio signify?

I did some research and according to Wikipedia an "enrichment factor" in terms of bioinformatics is a concept "to measure the added value of a search tool over another one or over the homogeneous distribution in the genome population." Does this mean that multitask DNN's are "better" than the other models, due to its higher number of EF/EFmax points?

But if we look at the curves on the right, it seems to say that for fewer number of compounds, the ratio for multitask DNN's is lower...

# Answer
> I did some research and according to Wikipedia an "enrichment factor" in terms of bioinformatics is...

Where possible, it's best to use explanations provided in the original source. This is because definitions could be different, depending on who is writing. Here's some descriptive information from the paper that may be more informative:

> The enrichment factor (EF) is another metric we use and is very popular for use in evaluation of virtual screening performance. EF denotes the quotient of true actives among a subset of predicted actives and the overall fraction of actives.... To evaluate the maximum achievable enrichment factor (EFmax), the maximum number of actives among the percentage of tested compounds is divided by the fraction of actives in the entire dataset. To more consistently quantify enrichment at 0.1% of all compounds tested, we report the ratio of EF/ EFmax. This normalized metric is useful because EF values are not directly comparable across different datasets due to the maximum possible EF being constrained by the ratio of total to active compounds and the evaluated fraction as shown in the equation above.

So for the purpose of this paper, they're saying that this ratio is an appropriately normalised substitute for a plain enrichment factor. High enrichment factor == good; low enrichment factor == bad.

Given the information you've provided, I'd agree that the curves in this figure appear to suggest that the multitask DNN performs worse with smaller numbers of active compounds (where Naive Bayes or Logistic regression would be a better choice). It's hard to tell from the graph (and legend) what the curves represent, and the clipping at EF/EFmax = 0 doesn't help.

The curves have a strange shape, suggesting the number of active compounds is dependent on the EF/EFmax ratio, rather than the other way round (which would make more sense to me). If I consider the curves to be something like a smoothed mean of the active compound number, then any place where there are fewer dots (i.e. low EF/EFmax for Multitask DNN) will have more associated error. For Multitask DNN, I see a large gap between 0.0 and 0.2, with lots of outliers on the 0.0 line (for all methods), so presumably the 0.0 outliers are skewing the curves.

Given that this is a preprint, the generation of the curves isn't specified, and the other figures (especially Figure 3) seem to tell a different story (i.e. that multi-task DNN performs better at all scales), I'd recommend contacting the authors to ask / clarify what the curves represent, or ignoring the curves entirely.
> 1 votes

# Answer
Ok I'll be brief, given the ms is likely under review.

**Background** in silico drug discovery has been part of deep learning for a long time 2015(?), and tensorflow was deployed. They worked out how to vectorise a small molecule and thereafter the entire field was.wide open to advanced ML.

**Ratio** The ratio is an accuracy index of the machine learning algorithms. The denominator is the total number of biological observations, the numerator is the predicted number. This 90% would be an accuracy index, i.e. the algorithm is correct 9 out of 10 times.

**Result** The plot shows, as the data set size increases the accuracy of the prediction increases most notably for the DNN and that conforms to all the theory. DNN/ANN only become powerful for 'big data', standard ML, i.e. random forest, naive Bayes, LR work well on small data sets, are fast and produce a classification which isn't wrapped up in a 3-4 Gig dump, thereafter performance saturates for routine ML methods.

The multi task DNN (MT-DNN) means that more than one loss function is being used, which is unusual in generic deep learning, but has been applied to natural language and drug discovery previously. Why this is being deployed is a long explanation, so I'll leave it there.
> 0 votes

---
Tags: machine-learning
---


# Question
Title: Changing this code in a way to work for my files
I have a bash script which extracts some information from a .vcf file but one .vcf file at each time. How I can change this script in a way to work on a bunch of .vcf files and the output is a .txt file named by the name of original .vcf

This is my script

help please

---

> EDITED

```
 #!/usr/bash

#outfilename
outname=$".parsed.txt"

#Header for output file
echo -e "$outname\t$chrom\t$Pos\t$Ref\t$Alt\t$ReadCount\t$VariantAlleleCount\t$ReferenceAlleleCount\t$ReadCountControl\t$VariantAlleleCountControl\t$ReferenceAlleleCountControl\t$VAF" >> $outname ;
while read -r line ;
do

#Basic information
chrom=$(echo "$line" | sed 's/ /\t/g' | cut -f 1) #&& echo $chrom;
Pos=$(echo "$line" | sed 's/ /\t/g' |  cut -f 2) #&& echo $Pos;
Ref=$(echo "$line" | sed 's/ /\t/g' | cut -f 4)
Alt=$(echo "$line" | sed 's/ /\t/g' | cut -f 5)

#Tumor sample read, variant and reference information
ReadCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 13 | sed 's/ReadCount=//' )
VariantAlleleCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 26| sed 's/VariantAlleleCount=//')
ReferenceAlleleCount=$(echo $ line | awk -v rc="$ReadCount" -v vac="$VariantAlleleCount" '{print rc-vac}')

#Control or Normal read, variant, reference information
ReadCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 14 | sed 's/ReadCountControl=//')
VariantAlleleCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 27 | sed 's/VariantAlleleCountControl=//')

ReferenceAlleleCountControl=$(echo "$line" | awk -v rcc="$ReadCountControl" -v vacc="$VariantAlleleCountControl" '{print rcc-vacc}')

VAF=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 28 | sed 's/VariantAlleleFrequency=//')
#Print output
echo -e "$outname\t""$chrom"'\t'"$Pos"'\t'"$Ref"'\t'"$Alt"'\t'"$ReadCount"'\t'"$VariantAlleleCount"'\t'"$ReferenceAlleleCount"'\t'"$ReadCountControl"'\t'"$VariantAlleleCountControl"'\t'"$ReferenceAlleleCountControl"'\t'"$VAF" >> $outname ;

#Remove info tags from VCF
done; < <( egrep -v '#' "$1")

```

---

And the error I obtained

```
[fi1d18@cyan01 new]$ bash vcf_parasing.sh
vcf_parasing.sh: line 7: syntax error near unexpected token `echo'
vcf_parasing.sh: line 7: `echo -e $outname'\t' "Chrom"'\t'"Position"'\t'"Ref"'\t'"Alt"'\t'"TumorReadCount"'\t'"TumorVariantAlleleCount"'\t'"TumorReferenceAlleleCount"'\t'"NormalReadCount"'\t'"NormalVariantAlleleCount"'\t'"NormalReferenceAlleleCount"'\t'"VAF" > $outname'
[fi1d18@cyan01 new]$

```

My script

https://www.dropbox.com/s/di49mij2a43i43s/vcf\_parasing.sh?dl=0

Two samples of .vcf

https://www.dropbox.com/s/ryk7q0aybeo5cy2/LP6008202-DNA\_B03\_vs\_LP6008201-DNA\_B03.snp.pass.vcf?dl=0

https://www.dropbox.com/s/u44s6d54r5d57bu/LP6008336-DNA\_H02\_vs\_LP6008333-DNA\_H02.snp.pass.vcf?dl=0

# Answer
> 5 votes
Using shell loops for text practice is considered bad practice. It is exceedingly slow, the syntax is complicated so it's very easy to get it wrong and it's just painful. The shell isn't designed as a proper scripting language, so while it can be (ab)used that way, you really should avoid it whenever possible.

Here's a Perl script that should do exactly what your bash script does but will be several times faster and can deal with multiple input files:

```
#!/usr/bin/env perl

## Iterate over all files given at the command line
foreach my $inFile (@ARGV) {
  print STDERR "Processing $inFile...\n";
  my $outFile = "$inFile.parsed.txt";
  $outFile =~ m|([^/]+)$|;
  #/#Useless comment to fix SE syntax highlighting
  my $outFileName = $1;
  open(my $outFileHandle, '>', $outFile) or
      die("Couldn't open '$outFile' for writing: $!");
  open(my $inFileHandle, '<', $inFile) or
      die("Couldn't open '$inFile' for writing: $!");

  ## Print header
  print $outFileHandle "$outFileName\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t" .
      "TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount" .
          "\tNormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n";
  ## Read all lines of the input file
  while (<$inFileHandle>) {
    ## Skip headers
    next if /^#/;
    ## Get VCF fields
    my ($chrom, $pos, $name, $ref, $alt, $qual,
        $filter, $info, $format, $values) = split(/\t/);
    my @infoFields = split(/;/, $info);
    my $readCount = $infoFields[12];
    $readCount =~ s/ReadCount=//;
    #/#Useless comment to fix SE syntax highlighting

    my $variantAlleleCount = $infoFields[25];
    $variantAlleleCount =~ s/VariantAlleleCount=//;
    #/#Useless comment to fix SE syntax highlighting
    my $referenceAlleleCount= $readCount - $variantAlleleCount;

    my $readCountControl = $infoFields[13];
    $readCountControl =~ s/ReadCountControl=//;
    #/#Useless comment to fix SE syntax highlighting

    my $variantAlleleCountControl = $infoFields[26];
    $variantAlleleCountControl =~ s/VariantAlleleCountControl=//;
    #/#Useless comment to fix SE syntax highlighting

    my $referenceAlleleCountControl = $readCountControl - $variantAlleleCountControl;
    my $VAF = $infoFields[27];
    $VAF =~ s/VariantAlleleFrequency=//;
    #/#Useless comment to fix SE syntax highlighting

    my $outString = join("\t", $outFileName, $chrom, $pos, $ref, $alt, $readCount,
                         $variantAlleleCount, $referenceAlleleCount, $readCountControl,
                         $variantAlleleCountControl, $referenceAlleleCountControl, $VAF);
    print $outFileHandle "$outString\n";
  }
  close($outFileHandle);
  close($inFileHandle);
}

```

Save that file as `parser.pl` and then run it on all your vcf files with:

```
perl parser.pl /path/to/vcfs/*vcf

```

That should create a `.parsed.txt` file for each input file in the target directory.

---

If you just have to use the shell for this, you can try with this improved version of your script:

```
#!/usr/bin/env bash

for file in "$@"; do
  outFile="${file}.parsed.txt"
  outFileName=$(basename "$outFile")
  (
    printf '%s\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t' "$outFileName"
    printf 'TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount\t'
    printf 'NormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n'
  ) > "$outFile"

  echo "Reading $file"
  grep -v '^#' "$file" |
    while read -r line; do
      read chrom Pos Ref Alt rest < <(awk -F"\t" '{print $1,$2,$4,$5}' <<<"$line")
      ReadCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 13 | sed 's/ReadCount=//' )
      VariantAlleleCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 26| sed 's/VariantAlleleCount=//')
      ReferenceAlleleCount=$(($ReadCount - $VariantAlleleCount))

      ReadCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 14 | sed 's/ReadCountControl=//')
      VariantAlleleCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 27 | sed 's/VariantAlleleCountControl=//')

      ReferenceAlleleCountControl=$(($ReadCountControl - $VariantAlleleCountControl))

      VAF=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 28 | sed 's/VariantAlleleFrequency=//')
      printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$outFileName" "$chrom" "$Pos" "$Ref" "$Alt" \
             "$ReadCount" "$VariantAlleleCount" "$ReferenceAlleleCount" "$ReadCountControl" "$VariantAlleleCountControl" \
             "$ReferenceAlleleCountControl" "$VAF" >> "$outFile"
    done
done

```

---

Just to illustrate what I mean about the shell being slow, I ran both scripts above on the example files you provided:

```
$ time foo.pl *vcf 2>/dev/null

real    0m0.051s
user    0m0.047s
sys     0m0.003s

$ time bar.sh *vcf 2>/dev/null

real    2m28.842s
user    3m33.648s
sys     1m12.456s

```

As you can see above, the perl script took less than a second, while the bash script took almost two and a half minutes!

---

Finally, if you find the perl script too long, here's the same thing in a more condensed style:

```
#!/usr/bin/env perl

## Iterate over all files given at the command line
foreach my $inFile (@ARGV) {
  print STDERR "Processing $inFile...\n";
  my $outFile = "$inFile.parsed.txt";
  $outFile =~ m|([^/]+)$|;
  #/#Useless comment to fix SE syntax highlighting
  my $outFileName = $1;
  open(my $outFileHandle, '>', $outFile);
  open(my $inFileHandle, '<', $inFile);

  ## Print header
  print $outFileHandle "$outFileName\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t" .
      "TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount" .
          "\tNormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n";
  ## Read all lines of the input file
  while (<$inFileHandle>) {
    ## Skip headers
    next if /^#/;
    ## Get VCF fields
    my ($chrom, $pos, $name, $ref, $alt, $qual,
        $filter, $info, $format, $values) = split(/\t/);
    $info=~s/[^;]+=//g;
    my @infoFields = split(/;/, $info);
    print $outFileHandle (join "\t", $outFileName, $chrom, $pos, $ref, $alt,
                               $infoFields[12],$infoFields[25],
                               $infoFields[12]-$infoFields[25],
                               $infoFields[13],$infoFields[26],
                               $infoFields[13]-$infoFields[26],
                               $infoFields[27]) . "\n";
  }
}

```

# Answer
> 1 votes
Here is a quick and kind of dirty solution with minimum intervention to the code:

* Add `$outname'\t'` right after `echo -e` to two two lines, the lines after `#Header for output file` and `#Print output`. The former will alter your header and the latter your actual output: you will have a first column with file name information.
* Then you can run this script in a for loop to generate multiple output files which then can be merged into a single file with `cat`.

Note that with this workaround, you will have multiple headers in the merged file. If you need a clean output, you will need to put your while loop in a for loop.

---
Tags: vcf, bash, linux
---


# Question
Title: How can I not show insertions in the Integrative Genome Viewer (IGV)?
I am using IGV 2.5.2 and would not like to see the insertions in my aligned reads. How can I do that? I have managed to remove mismatched bases because there is an option to do so when I right-click on the window ("Show mismatched bases"). Now I would like to do the same with these purple insertions.

### 3rd January 2020 update

Thanks to a member's response, I found some parameter in the preferences that might be helpful to hide the indels. However, I have tried different numbers from 0 to 1000 and the indels are still shown. What numbers should I put to the parameters "Show indel threshold (bases)" and should I click the box or not?

# Answer
Preferences `->` Alignments `->` Hide indels \< 'x' bases.

---

Credit to Pierre Lindenbaum
> 3 votes

---
Tags: sequence-alignment, visualization, igv
---


# Question
Title: How to find which scales are available in SeqUtil of Biopython?
I am trying to analyse the hydrophobicity of a sequence using BioPython's "SeqUtil".

```
analysed_seq = ProteinAnalysis("AKVLGSLEALEPVHYEEKNWCE")
analysed_seq.protein_scale(Scale, WindowSize, Edge)

```

**Known parameters**

* `WindowSize` refers to the total number of amino acids to be considered either side of the residue (an `int` variable).
* `Edge` refers to the weighting of the scale further from the centre of the window (0-1). \`

**Unknown parameter**

* Scale`refers to which scale being used. I am interested in how to give the`Scale\` variable a valid variable.

The Biopython docs simply point out that there are many `scales` to choose from and the source code doesn't seem to give any clues about scale options.

Where can I find a list of the `scale` are available and how to call them?

# Answer
> 1 votes
Use this line of code:

```
print(analysed_seq.protein_scale(analysed_seq.get_amino_acids_percent(), window=9, edge=0.4))

```

---
Tags: phylogenetics, python, phylogeny
---


# Question
Title: mapping and de novo assembly of plant without reference genome
I'm studying on a plant that has no reference genome and only has one scaffold assembly and one gff3 annotation file. Can I create an index with the same assembly and gff3 in STAR and do the mapping? If I want to do de novo assembly using the Galaxy Can I import the links of 20 fastq files from the RNA-seq of the samples together to construct the Reference Genome with Trinity? If this is possible, please help me import the links of these files together into Galaxy and make the reference genome.

Thanks a lot

# Answer
You have multiple questions, I will try to answer one by one:

* On creating an index with `STAR`: As long as you provide a fasta file and its corresponding gff, you will be able to generate an index.
* On creating a reference with `Trinity`: You have mentioned "to construct the Reference Genome with Trinity", however, with your RNA-seq data and Trinity, you can construct a de novo transcriptome and not a genome.
* Uploading files to `Galaxy`: I suggest to do so via FTP using `FileZilla` or something similar. The transfer can take some time if your files are large though.

If I understand your questions correctly, you are performing RNA-seq and you are concerned about the completeness of the genome of your model organism. In that case, I would advise to bring together a de novo transcriptome from your RNA-seq data (with Trinity), use one of the quasi-mapping methods (with a tool like `Salmon` or `kallisto`).
> 1 votes

---
Tags: rna-seq, assembly
---


# Question
Title: bedtool intersect to compare .vcf file to reference genome
I have a `.vcf` file generated using `samtools`. I want to know mutations in specific genes by comparing `.vcf` files with a reference genome. I cannot use the Ensembl VEP tool as my reference genome is not listed there. Can you please help?

# Answer
> 1 votes
If you want to filter a VCF file, you might want to give **VCFtools** a try.

Here is the documentation : VCFtools Manual

I am not a specialist of this tool but you could do as the following. First build a bedfile with the chromosome names and positions of your genes (start and stop) such as this :

```
chr   start   stop   (name)   (score)   (strand)

```

Then it should work with this command line :

```
vcftools --vcf file.vcf --bed positions.bed --out filtered.vcf

```

The output will be a VCF filtered on the provided regions

---
Tags: vcf
---


title: Looping over several files in bash

text:
# Question
Title: Looping over several files in bash
I have a bash script:

I am wondering how I can change this script to loop over a bunch of .vcf files and give output .txt with the name of corresponding .vcf

I tried changes done in similar script in Changing this code in a way to work for my files but I just messed everything up

# Answer
> 2 votes
If we save the script pasted in the main post as a sh file and we have some .vcf files in a folder, by this line we can iterate over vcf files to extract what mentioned in the script returning .txt output by the name of corresponding vcf

`for file in *.vcf ; do bash indel_vcf_parasing.sh $file ; done`

# Answer
> 1 votes
Error means that quotations are missing. I am not going to review this chunk but at first glance check

```
printf 'TumorAltAlleleCount\tTumorAltAlleleFrequency\t'\n'

```

There are three quotation marks so one is not matched. I suggest you develop code within a good IDE that will help you debug. Please google around for potential solutions. I personally do R and Shell scripting in RStudio. It differently colors quoted and unquoted code so you immediately see if you missed a closing quotation somewhere.

---
Tags: vcf, bash, linux, wgs
---


# Question
Title: Single-cell sequencing dataset has too many barcodes
I am analyzing a single-cell sequencing dataset from the website 10xgenomics, with 2000 cells. It is a BAM file and I am trying to obtain the individual cells per sample. I used the command

```
samtools view mtdnaAsorted.bam | grep CB:Z: | sed 's/.CB:Z:([ACGT]).*/\1/' | sort | uniq -c > reads_per_barcode_mtdnaA

```

to get a list of the individual bar codes and their reads, but there are 100,000 unique barcodes even with only 2000 cells. Does anyone know what could be the problem?

I know that there could be noise, but this is a problem for me because I'm trying to analyze coverage per cell, and if there are a lot of cells that contribute towards the total number of reads, the net coverage per cell is going to be really low

# Answer
> 4 votes
The BAM file you downloaded has all of the detected 10X barcodes instead of those that represent cells. In a given 10X experiment the number of input barcodes vastly outnumbers the input cell count. After the reaction some of barcoded beads are invariably sequenced because they are encapsulated with ambient RNA, dead cells etc. After processing the BAM file Cell Ranger infers which cells likely represent cells and generates gene expression matrices from these files.

You would need to download one of the output files that contains the list of passed barcodes such as the `Gene / cell matrix (filtered)` and then use the list of passed barocode sequences to filter your BAM file. Alternatively you could approximately select the same cells by selecting the top 2,000 barcodes with the highest numbers of reads.

---
Tags: samtools, single-cell, 10x-genomics
---


title: Why don't all reads have adaptors?

text:
# Question
Title: Why don't all reads have adaptors?
I got NGS reads back from sequencing platform. I check for adaptors and trimmed them. But I realized only a fraction (eg 30%) have adaptors... why not all of them?

Thanks

# Answer
Things like this might depend on your specific library prep, but in general: sequencing starts at the end of the adapter, not before it. You will only see adapters if you sequence through the entire fragment into the adapters on the other side. If your fragment is long enough you won't see any adapter, and of course, that's the desirable outcome of your library prep.
> 8 votes

---
Tags: ngs, adapter, cutadapt, fastqc
---


# Question
Title: Downloading multiple SRA files from several SRA accession IDs does not work
I am trying to download multiple SRA files located in several SRA accessions. Some of my accession numbers are as follows:

```
SRR1002678
SRR1158898
SRR1158907
SRR1158950
SRR1181152
SRR1181155
SRR1181196
SRR1181204
SRR1181205
SRR1181208
SRR1269497
SRR2100227
SRR2100419
SRR2100830
SRR2101344
SRR671725
SRR671726
SRR671740
SRR671759
SRR671760
SRR671761
SRR671779
SRR671782
SRR671815
SRR671821
SRR671842
SRR671862
SRR671863
SRR671864
SRR671865
.............

```

I tried the following command:

```
for i in {SRR1002678..SRR1184123}; do wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR1*/$i; done

```

but this did not work.

I tried `prefetch` as well:

```
sort -u SRR.txt | parallel -j 4 "prefetch {}"

```

where SRR.txt is the list of my SRA IDs. However, with this command only a few files are downloaded as \*.sra files, while with several other SRA IDs I get the following error:

```
2018-11-15T06:17:20 prefetch.2.8.2: 1) Downloading 'SRR1183962'...
2018-11-15T06:17:20 prefetch.2.8.2:  Downloading via https...
2018-11-15T06:20:21 prefetch.2.8.2: 1) 'SRR1183962' was downloaded successfully
2018-11-15T06:20:27 prefetch.2.8.2: 'SRR1183962' has 0 unresolved dependencies
2018-11-15T06:20:38 prefetch.2.8.2: 'SRR1183962' has no remote vdbcache

2018-11-15T06:18:59 prefetch.2.8.2: 1) Downloading 'SRR1184034'...
2018-11-15T06:18:59 prefetch.2.8.2:  Downloading via https...
2018-11-15T06:20:40 prefetch.2.8.2: 1) failed to download SRR1184034
2018-11-15T06:20:40 prefetch.2.8.2 sys: libs/kns/unix/syssock.c:606:KSocketTimedRead: timeout exhausted while reading file within network system module - mbedtls_ssl_read returned -76 ( NET - Reading information from the socket failed )
2018-11-15T06:20:40 prefetch.2.8.2 int: libs/kns/unix/syssock.c:606:KSocketTimedRead: timeout exhausted while reading file within network system module -

```

I do not have a range of SRA files to download but quite random SRA files to download. What am I doing wrong? Any help is highly appreciated.

# Answer
> 6 votes
First of all, the brace expansion you tried to use isn't valid. While bash can expand `{1..3}` to `1 2 3`, it has no way of knowing how to expand `{SRR1002678..SRR1184123}`:

```
$ echo {1..3}
1 2 3
$ echo {SRR1002678..SRR1184123}
{SRR1002678..SRR1184123}

```

Perhaps you meant `SRR{1002678..1184123}`? But there's no need for this, you already have the accessions in a file and can read that. I see that the FTP server you are accessing has organized the data like this:

```
sra/sra-instant/reads/ByRun/sra/SRR/SRRNNN/SRRNNNXXXX/SRRNNNXXXX.sra

```

Where `NNN` are the 1st three digits of the accession. So, for example, the file for accession `SRR1002678` is at:

```
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR118/SRR1181152/SRR1181152.sra

```

With this in mind, you can write a little script that extracts the URL from each accession:

```
while read i; do
    n=$(echo "$i" | sed -E 's/(SRR...).*/\1/');
    wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/$n/$i/$i.sra
done < SRR.txt

```

The only trick is the `sed` command which simply looks for `SRR` and the three next characters and removes everything else on the line (it "captures" the found pattern and then replaces the line with what was captured (`\1`)).

# Answer
> 0 votes
For me it works best to combine `cat`, `while` and official download tool `fasterq-dump`.

```
cat your_accession_list.txt | while read i
do
fasterq-dump ${i} <options>
done

```

As commented above, `your_accession_list.txt` should have Linux endings. Link to `fasterq-dump` is here.

---
Tags: linux, sratoolkit
---


# Question
Title: Infer the new (Ban et. al) ribosomal nomenlature (ex. uL53 ) from the ribosomal protein's sequence
##

I'm trying to programmatically construct a name for each protein subchain in any ribosome from Uniprot in accordance with Ban et. al's 2014 proposal (excerpt given below) using PDB's and Uniprot's APIs. This is to be a part of a bigger ribosomal analysis suite released next year and we would like to standardize the names as much as possible. Ban et. al's proposal: \[https://bangroup.ethz.ch/research/nomenclature-of-ribosomal-proteins.html\]

---

I'm however not remotely a biologist (I work in compilers) and ***would really appreciate some clarification on the way nomenclature works across domains of life and how homology plays into it.*** This is my *naive* approach thus far:

* I get a ribosome from PDB let's say (*3J9M* or *5MYJ* or any other for that matter) and split it into protein subchains. So on the order of 20-80 proteins.
* For each protein subchain, I can get all the information in the world on it from Uniprot: names, sequence, organisms containing it, publications, synonymous names, you name it...
* Right now (again, naively) for those proteins that *don't already have a new name assigned to them*, I take name-synonyms that PDB has plenty of (from various publications, I assume) and scan them with regexes (ex.`/[LS]\d{1,2}/g`) for things like `"L8"` or `"S15"`.
* Thus, for each subchain, I end up with something like `"L15"` or `"S28"` in multiple copies (depending on the abundance of synonymous entries in PDB).
* The rut i'm in right now is this: having something like `"L15"` to match it to a new "Ban-nomenclature" name: they provide the nomenclature look-up tables in terms of these three "historical" naming conventions(i.e "human", "yeast", "bacteria") and a taxonomic range while all the I have is the ribosomal structure(ex \`\` from which I got the subchain. I'm sure there's plenty of information that I can get

Could somebody elucidate to me which assumptions I can and cannot make in this field and how exactly to use the *taxonomic range*?

Say, (i) can I assume that all archean names are the same as bacterial homologs? If I draw a eukaryotic cytoplasmic ribosome, do I match its subchains against yeast column or human column (since both are eukaryotic)?

---

Somebody suggested that I start grouping these subchains based on sequence, not on a domain and I do have access to each protein's sequence, but I am not sure how to go about it. I was thinking to get each of Ban's new names' **UniRef90** cluster, let's say, and then start checking which cluster a given subchain belongs to.

But I'm obviously guessing at this point. Any advice or pointers are very appreciated.

##

# Answer
*This is a copy-paste nearly verbatim of comments so that the question has an answer*

What you are trying to do is find what genes are in what homology cluster. This is common problem and there are many solution each with some issues.

Uniprot90 is indeed a cluster of homologues but it is too limited. Whereas you require clusters that span all the domains (universal).

### COG

The nearly original grouping, is COG (cluster of orthologous genes). It is really handy for this but aren't really available from Uniprot and NCBI (its curator) has been phasing COGs out. In Uniprot records, they are *patchily* present as its derivative eggNOG (which has varying degrees of "scope", LUCA being the universal one):

```
<dbReference type="eggNOG" id="COG0451">
<property type="taxonomic scope" value="LUCA"/>
</dbReference>

```

### Pfam

Another option that is well mapped and always present when possible in Uniprot entries is Pfam ids. These are domain family annotations and are also good option. It is a *domain* fold group. So domains within protein with homologous folds are grouped into one. If a protein has two known domains, it will have two Pfam entries —repeat proteins have loads, but most protein have one domain. In Uniprot XML you will find these are:

```
<dbReference type="Pfam" id="PF01370">

```

Similarly there is InterPro, which is more narrow than Pfam, but generally overlaps.

### DIY

Then there's always the option of making your own psi-Blast motif on the NCBI website (setting high number of hits, setting the database to refseq or even PDB and doing a few iterations), saving the motif pattern (PSSM) and using local blast with it.

### Caveats

In general you may however encounter these pitfalls:

1. you may have two protein clusters that share a single id
2. you have gene fusions (one protein, 2 domains of interest). In Pfam the protein will have two ids, with COGs or similar, you'd have only an anamalous length to go by.
3. you may have wished to find analogues (i.e. two protein that arose independently that have the same role (a lot of ribosomal protein are simply structural). No homology clustering scheme can find these. Say a bacterium has evolved to have a different protein with that role, yet as it has never been studied it is a domain of unknown function (guilt by association is the common strategy, but that is beyond your scope here).
4. you have a very small protein which is poorly conserved
> 1 votes

---
Tags: python, uniprot, ribosomal, nomenclature, sequence
---


title: Merging matched parts of two dataframes

text:
# Question
Title: Merging matched parts of two dataframes
I have a big data frame like this

```
> dim(data1)
[1] 474303    128
> head(data1[1:2,1:6])
   Tumor_Sample_Barcode Chromosome Position
1:                   A1       chr1       15772155
2:                   A1       chr1       16054813
   End_Position Reference_Allele Tumor_Seq_Allele2
1:     15772155                G                 A
2:     16054813                G                 A
>

```

I have another dataframe like this

```
Chrom   Position    Ref Alt VAF
chr1    2716062 C   A   0.323
chr1    3095112 T   G   0.365
chr1    3342771 C   G   0.123
chr1    3526376 C   T   0.365
chr1    4121415 T   G   0.225
chr1    4215706 G   A   0.125
chr1    4304014 T   G   0.286
chr1    4304016 A   C   0.302
chr1    4388809 A   C   0.135

```

I want to create another column named ***i\_TumorVAF\_WU*** in my big data frame in which VAF column comes from small data frame. Please notice that VAF values should be matched with Position in big data because small data coming from big data

Can you help please?

```
    > unique(merge(data1,data,x.by=data1$Position,y.by=data$Position))

```

Returns a weird output with losing a big part of big data frame. I just want to add corresponding VAF values to the big data by keeping all the information in big data file

# Answer
> 3 votes
If you have issue with memory and dealing with large object, maybe `data.table` is the way to go (https://github.com/Rdatatable/data.table/wiki):

Let's take two dummy data.frame:

```
dfA <- data.frame(Chrom = rep("chr1",150),
                  Pos = 1:150,
                  ref = rep("A",150))
dfB <- data.frame(Chrom = rep("chr1",10),
                  Pos = 11:20,
                  Vaf = rep("VAF", 10))

```

The big data.frame `dfA`:

```
> head(dfA)
  Chrom Pos ref
1  chr1   1   A
2  chr1   2   A
3  chr1   3   A
4  chr1   4   A
5  chr1   5   A
6  chr1   6   A

```

and the small dataset:

```
> dfB
   Chrom Pos Vaf
1   chr1  11 VAF
2   chr1  12 VAF
3   chr1  13 VAF
4   chr1  14 VAF
5   chr1  15 VAF
6   chr1  16 VAF
7   chr1  17 VAF
8   chr1  18 VAF
9   chr1  19 VAF
10  chr1  20 VAF

```

Now, using `data.table` that is more efficient and less memory consuming for the manipulation of dataframe, you can do:

```
library(data.table)
setDT(dfA)
setDT(dfB)
dfA[dfB, on = c("Chrom","Pos"), vaf := i.Vaf]

```

And now, your data.frame `dfA` looks like:

```
> dfA[5:25]
    Chrom Pos ref  vaf
 1:  chr1   5   A <NA>
 2:  chr1   6   A <NA>
 3:  chr1   7   A <NA>
 4:  chr1   8   A <NA>
 5:  chr1   9   A <NA>
 6:  chr1  10   A <NA>
 7:  chr1  11   A  VAF
 8:  chr1  12   A  VAF
 9:  chr1  13   A  VAF
10:  chr1  14   A  VAF
11:  chr1  15   A  VAF
12:  chr1  16   A  VAF
13:  chr1  17   A  VAF
14:  chr1  18   A  VAF
15:  chr1  19   A  VAF
16:  chr1  20   A  VAF
17:  chr1  21   A <NA>
18:  chr1  22   A <NA>
19:  chr1  23   A <NA>
20:  chr1  24   A <NA>
21:  chr1  25   A <NA>
    Chrom Pos ref  vaf

```

As you see, we incorporate the column `Vaf` from the small dataframe into the big dataframe.

With your example, I will first rename the colnames of your small dataframe in order to have the perfect match with the big dataframe:

```
colnames(data) <- gsub("Chrom", "Chromosome", colnames(data))

```

then,

```
setDT(data1)
setDT(data)
data1[data, on = c("Chromosome","Position"), i_TumorVAF_WU := i.VAF]

```

Does it look what you are trying to obtain ?

# Answer
> 1 votes
Read the manuals for `data.table::merge` and `base::merge`.

This (wrong syntax):

```
unique(merge(data1,data, x.by = data1$Position, y.by = data$Position))

```

Should work fine as below:

```
merge(data1, data, by = "Position")

```

Regarding:

> "losing a big part of big data frame"\`

It is only returning rows when there is a match on *Position*, if you need to return the full dataset, then use:

```
merge(data1, data, by = "Position", all.x = TRUE)

```

---
Tags: r, wgs, maf, snv
---


# Question
Title: Changing the name of elements in a column
I have a very big matrix like

```
> head(mat[1:4,1:4])
                     Tumor_Sample_Barcode Chromosome Start_Position End_Position
1: LP6005334-DNA_H01_vs_LP6005333-DNA_H01       chr1        1026918      1026918
2: LP6005334-DNA_H01_vs_LP6005333-DNA_H01       chr1       15772155     15772155
3: LP6005334-DNA_H01_vs_LP6005333-DNA_H01       chr1       16054813     16054813
4: LP6005334-DNA_H01_vs_LP6005333-DNA_H01       chr1       24120222     24120222

> dim(mat)
[1] 450754     53
>

```

In `Tumor_Sample_Barcode` column I have the name of samples so that

```
> unique(mat$Tumor_Sample_Barcode)
 [1] "LP6005334-DNA_H01_vs_LP6005333-DNA_H01"     "LP6005500-DNA_D03_vs_LP6005499-DNA_D03"
 [3] "LP6005935-DNA_G02_vs_LP6005979-DNA_C01"     "LP6007600_vs_LP6007599"
 [5] "LP6008202-DNA_B03_vs_LP6008201-DNA_B03"     "LP6008269-DNA_F06_vs_LP6008268-DNA_G06"
 [7] "LP6008334-DNA_A03_vs_LP6008335-DNA_E02"     "LP6008334-DNA_A04_vs_LP6008335-DNA_E03"
 [9] "LP6008334-DNA_B02_vs_LP6008335-DNA_F01"     "LP6008334-DNA_C02_vs_LP6008335-DNA_G01"
[11] "LP6008334-DNA_D02_vs_LP6008335-DNA_H01"     "LP6008336-DNA_F02_vs_LP6008333-DNA_F02"
[13] "LP6008336-DNA_G01_vs_LP6008333-DNA_G01"     "LP6008336-DNA_H01_vs_LP6008333-DNA_H01"
[15] "LP6008336-DNA_H02_vs_LP6008333-DNA_H02"     "LP6008337-DNA_A07_vs_LP6008338-DNA_E06"
[17] "LP6008337-DNA_H06_vs_LP6008338-DNA_D06"     "LP6008460-DNA_A04_vs_LP6008340-DNA_D05__pv"
[19] "LP6008460-DNA_D01_vs_LP6008340-DNA_E02__pv" "LP6008460-DNA_F02_vs_LP6008340-DNA_E05"
[21] "LP6008460-DNA_G03_vs_LP6008340-DNA_C05__pv"
>

```

I want to change the name of these samples to A1... A21 but I don't know how

Can you help?

Thank you

# Answer
> 3 votes
If your data is in a dataframe you can make use of `plyr`. Using `mapvalues()` from this library, you can rename provide a list of keys and a list of values. Instead of explicitly listing the keys, we can just use `unique(mat$Tumor_Sample_Barcode)` \[credit to @dc37\]. I use the `paste0` function here and the colon operator to generate the vector of `A1,A2,...,A21`.

Assign the vector produced by `mapvalues()` to the column and you're off!

```
install.packages("plyr")
mat$Tumor_Sample_Barcode <-
plyr::mapvalues(x=mat$Tumor_Sample_Barcode,
                from=unique(mat$Tumor_Sample_Barcode),
                to=paste0("A", 1:21))

```

I'd recommend keeping this column and adding a new column with your simplified sample names, but this is not what you requested.

# Answer
> 2 votes
Change to *factor* then to *numeric*:

```
mat$Tumor_Sample_Barcode <- as.numeric(as.factor(mat$Tumor_Sample_Barcode))

```

Using @dc37 example data:

```
as.numeric(as.factor(df$vect))
# [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

```

If we need to prefix with "A"

```
paste0("A", as.numeric(as.factor(df$vect)))
# [1] "A1" "A2" "A3" "A4" "A1" "A2" "A3" "A4" "A1" "A2" "A3" "A4" "A1" "A2" "A3" "A4" "A1" "A2" "A3" "A4"

```

Or we can keep them as *factors*:

```
x <- sort(unique(df$vect))
factor(df$vect, levels = x, labels = paste0("A", seq_along(x)))
# [1] A1 A2 A3 A4 A1 A2 A3 A4 A1 A2 A3 A4 A1 A2 A3 A4 A1 A2 A3 A4
# Levels: A1 A2 A3 A4

```

# Answer
> 1 votes
Using `sub` in `base` R, you can:

1) define the replacement vector:

```
vector = paste0("A",seq(1:21))

```

2) and loop over each levels of your first columns to replace them individually:

```
for(i in 1:length(levels(mat$Tumor_Sample_Barcode)))
{
  levels(mat$Tumor_Sample_Barcode)[i] = sub(levels(mat$Tumor_Sample_Barcode)[i],vector[i], levels(mat$Tumor_Sample_Barcode)[i])
}

```

**Example:**

Here, a small example to show you:

```
df = data.frame(vect = factor(rep(c("a","b","c","d"),5)),
                       value = rnorm(20))

> head(df)
  vect        value
1    a  0.418186287
2    b -0.232936122
3    c  1.563349091
4    d -0.001453912
5    a -0.453995163
6    b -0.481360070

```

```
vector = paste0("A",seq(1,4))
for(i in 1:length(levels(df$vect)))
{
  levels(df$vect)[i] = sub(levels(df$vect)[i],vector[i], levels(df$vect)[i])
}

```

And you get:

```
> head(df)
  vect        value
1   A1  0.418186287
2   A2 -0.232936122
3   A3  1.563349091
4   A4 -0.001453912
5   A1 -0.453995163
6   A2 -0.481360070

```

Does it answer your question ?

---
Tags: r
---