# Question
Title: pymol script to only select 5 chains within a distance from a reference and save the selection
I am generating 30 asymmetric units of a protein structure. Now, I want to select different combinations of 5 units which are 30 Angstroms away from the reference and save it as a different file. Here is the command I could work out till now :

```
fetch 3q2v
symexp sym,3q2v,(3q2v),30
select chains, bychain all within 30 of 3q2v

```

How to loop in and select only 5 chains at a time obtained from the setup obtained by the above code, and save it in separate files?

# Answer
> 1 votes
Firstly, I would preface that the PyMOL command line in the GUI is great, but **you can run PyMOL as a regular Python library** (calling the appropriate python). This makes life so much simpler for more complex manipulations. You can do really complex manipulations with it. I cannot suggest it enough. The API version of the commands works only marginally differently, but you can do API calls in the GUI `cmd.`) and PyMOLWiki documentation is only slightly weaker for them (`help()` fills the gaps).

But if you really want to use the GUI route...

Let's get the names of the models in the selection

```
onames = cmd.get_object_list('chains')
print onames

```

Now, you want to save subsets of 5 chains (technically, objects or models). That is, you do not care about the order they are in and you don't want the same one twice (actually selections with 4 or less), so you want all combinations. Make sure to do the following print len(onames)

And google that number as `<n> choose 5`. That is how many files you will get. I am pretty sure most of these sets will be meaningful and will be useless. In a script a loop is easy, while in the gui you have to keep it simple or use a list coprehension hackishly.

```
 from itertools import combinations
 for s in combinations(onames, 5): cmd.save('-'.join(s)+'.pdb', ' or '.join(s))

```

Obviously, you could alternatively use an enumerate method around the iterator and get an index for the file. But do note f-strings will not work.

```
 for i, s in enumerate(combinations(onames, 5)): cmd.save(str(i)+'.pdb', ' or '.join(s))

```

---
Tags: protein-structure, visualization, pymol
---


# Question
Title: Common QC parameters and threshold in human WGS pipeline
I am trying to summarize a list of QC-parameters and its threshold that should be applied for human WGS in general to avoid bad quality of data.

Which qc parameters for WGS are you using and why?

Thanks!

# Answer
> 2 votes
I would suggest you look at the website for the excellent MultiQC tool (https://multiqc.info/). We use it in a clinical setting to collate all the QC statistics from our pipelines into a single report (It displays results from multiple QC programs such as FastQC, Picard, bcl2fastq, VerifyBAMid etc). The homepage has a selection of example reports - clicking on "Whole Genome Seq" shows a report with the metrics most labs would expect from their WGS pipeline.

If you dig around in their GitHub repo (https://multiqc.info/) you will get a good idea of the default thresholds set for flagging results - there should be a config file somewhere specifying all these values.

# Answer
> 0 votes
For your negative controls, you might want to check out QC Fail (https://sequencing.qcfail.com/). There you can see examples of how things look in QC when things have gone wrong, such as loss of quality, position-specific errors on flowcells and contamination with primer dimers.

---
Tags: wgs, qc
---


# Question
Title: How can I export a full alignment from IGV as an image?
The IGV browser lets you export an alignment as an image (File =\> Save Image). However, this image only contains those reads that fit in the viewing window:

As you can see in the image above by looking at the scroll bar on the right, this is only a small subset of aligned reads. Is there any way to export an image of the aligned reads that includes *all* reads and not only those currently shown in the window?

# Answer
> 4 votes
Running an IGV batch script helps increase the number of reads per image, compared with `File | Save image`, but it does not solve the issue by itself.

To run a script, first make a script file like so with 1 command per line:

```
snapshotDirectory ./
goto chr2:28582094-28582790
maxPanelHeight 10000
snapshot chr2_28582094_28582790.png

```

Run the script using `Tools | Run batch script` menu. The importat trick is the `maxPanelHeight 10000` which is what makes this include more reads than the default snapshot.

You can also convert this into a little shell script:

```
#!/bin/bash
mkdir -p /tmp/igv
sleep 2
echo "snapshotDirectory /tmp/igv"
sleep 2
goto chr2:28582094-28582790
sleep 2
echo "maxPanelHeight 10000"
echo "snapshot chr2_28582094_28582790.png"

```

Then, just open IGV, load the relevant BAM file and run the script like this to connect to IGV via telnet (credit to Pierre Lindenbaum):

```
sh foo.sh |  telnet 127.0.0.1 60151

```

---
Tags: sam, ngs, visualization, igv
---


title: SRA Toolkit execution problem

text:
# Question
Title: SRA Toolkit execution problem
I am trying to install SRA Toolkit on MacOS, following all the instructions from NCBI website. It is installed but fastq-dump cannot be executed because "cannot execute binary files". I tried "chmod +x", noah from home-brew to make binary executable on MacOS but fast-dump still cannot be executed. I attached the screenshot of my terminal window. Can anyone help?

# Answer
> 1 votes
You can't run Linux binaries on OSX (or vice versa). Download Mac OSX binaries.

---
Tags: fastq, sratoolkit
---


title: Pattern mining from a genomic sequence

text:
# Question
Title: Pattern mining from a genomic sequence
I need to find the following pattern from a genomic sequence

```
TTCGB-N15-TTCGB
B represents C, G, or T
N15 represents A, T, G or C and these four base combination should have the length of 15 bases.

```

To do the same, I tried `grep` command, but I could not specify `B` and `N15` in the `grep` command line. Therefore, please help me to find the same.

For example, I have a nucleotide sequences as follows

```
>ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT

```

Expected outcome,

> TTCGTGCATGCATATGCGCATTCGC

Thanks in advance.

# Answer
```
echo ">ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT" | grep -oP  'TTCG[CGT][ATGC]{1,15}TTCG[CGT]'

```
* `[CGT]` looks for one occurrence of one of the three.
* `[ATGC]{1,15}` looks for "up to 15 letter long" combinations of the four bases.

Moreover, `-P` makes your regex call `Perl-like` and is required to make this regex work. Different regex flavors behave (slightly) differently.
> 4 votes

---
Tags: r, python, bash, perl, shell
---


# Question
Title: How do researchers account for heterozygosity in genome assembly?
this is a very general question about diploid genome assembly. I am wondering how people deal with heterozygous region in a genome when assembling ? Do they pick one of the haplotype to be in the referene and leave out the other ? I am particularly thinking about situation similar to mammalian chromosome Y. Do you put Y or X in the assembly ? An hybrid ?

Thanks

# Answer
Generally a single allele is used. This may be arbitrary or match the major allele frequency of the species. For sex chromosomes, usually both/all are included, so for mammals both X and Y are part of the reference genome.
> 1 votes

---
Tags: reference-genome, genome-sequencing
---


title: Replacing some strings in R or terminal

text:
# Question
Title: Replacing some strings in R or terminal
I have a .txt file; I have 139 unique sample IDs in sampleID column

```
> head(mutations)
                                                                              sampleID chr     pos ref mut
4 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1  920552   G   A
5 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1350208   G   A
6 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1447367   T   G
7 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1909310   G   A
8 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 2172675   G   C
9 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 2846623   C   T
>

> dim(mutations)
[1] 3574407       5
>

```

I want to change name of sample from 1 to 139 in R or Linux terminal but I don't know how to do that, any help please?

Thanks

# Answer
> 1 votes
I believe that

```
uniqsampleID = 1:length(unique(mutations$sampleID))
names(uniqsampleID) = unique(mutations$sampleID)
mutations$sampleID = uniqsampleID[mutations$sampleID]

```

should work, assuming `mutations` is a data frame.

---

Explanation:

`1:N` creates a vector of length N, containing 1, 2, ..., N. So using `1:length(X)` will do this for each item in X, so `1:length(unique(mutations$sampleID))` should give you a vector of `1:139`.

Then `names(uniqsampleID) = unique(mutations$sampleID)` assigns to each of these numbers one of these sample IDs, which allows you to retrieve the sample number with `uniqsampleID[mutations$sampleID]`.

---
Tags: r, linux, looper
---


title: How I loop over this code

text:
# Question
Title: How I loop over this code
I have 139 samples; In R for each sample I run below codes

for example for ***sample 1***

```
sample_1 = whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = 1,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')

a1=(sample_1[[1]])

> class(sample_1)
[1] "list"

> a1=(sample_1[[1]])

> head(a1)
  Signature.1A Signature.1B Signature.2 Signature.3 Signature.4 Signature.5 Signature.6 Signature.7 Signature.8
1    0.5370653            0           0           0           0           0           0           0           0
  Signature.9 Signature.10 Signature.11 Signature.12 Signature.13 Signature.14 Signature.15 Signature.16 Signature.17
1           0            0            0            0            0            0            0            0    0.3383947
  Signature.18 Signature.19 Signature.20 Signature.21 Signature.R1 Signature.R2 Signature.R3 Signature.U1 Signature.U2
1    0.1123409            0            0            0            0            0            0            0            0
> dim(a1)
[1]  1 27

```

and for ***sample 2***

```
sample_2 = whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = 2,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')

a2=(sample_2[[1]])

```

And finally binding `a1,a2,...a139` by

```
sig=rbind(a1,a2)

> print(sig)
  Signature.1A Signature.1B Signature.2 Signature.3 Signature.4 Signature.5 Signature.6 Signature.7 Signature.8
1    0.5370653            0           0           0           0           0           0           0           0
2    0.7162292            0           0           0           0           0           0           0           0
  Signature.9 Signature.10 Signature.11 Signature.12 Signature.13 Signature.14 Signature.15 Signature.16 Signature.17
1           0            0            0            0            0            0            0            0    0.3383947
2           0            0            0            0            0            0            0            0    0.2837708
  Signature.18 Signature.19 Signature.20 Signature.21 Signature.R1 Signature.R2 Signature.R3 Signature.U1 Signature.U2
1    0.1123409            0            0            0            0            0            0            0            0
2    0.0000000            0            0            0            0            0            0            0            0
>

```

`sigs.input` is a common R object here

Doing these manually is too painful; Can you please help me how I do that automatically in R? I am not sure but something like looping or whatever you guess possible

Thanks

# Answer
```
list.of.results <- lapply(seq(INSERT_TOTAL_NUMER_OF_SAMPLES_HERE),function(sample.id){
whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = sample.id,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')[[1]]
})

sig <- do.call(rbind,list.of.results)

```
> 6 votes

---
Tags: r, loop
---


title: Method to detect genome doubling

text:
# Question
Title: Method to detect genome doubling
CNV is usually detected using either read depth or SNP VAF. However, neither of these methods can detect CNV in case of genome doubling (assuming genome doubling is a kind of genome scale CNV comparing to normal tissue which is diplod). We thought of other platform such as SNP array and we did not think they can detect genome doubling as well. Can any one share some comments on the method for such event? Thanks in advance.

# Answer
> 3 votes
In a tissue or cell line that used to be diploid and has undergone genome wide doubling, you expect more variants having observed allele frequencies that would also occur in tetraploid organism. That is, observed VAF distribution will have (possibly small) peaks at 0.25 and 0.75, in addition to the diploid peaks at 0, 0.5 and 1. Compare the distributions in the tissue or cell line of interest to those of the control, that is, known diploid. Also if it is a doubled male genome, the sex chromosomes will show a new peak at AF 0.5, that is, heterozygous X or Y variants. This is just an idea, I have never tried to look for such events.

# Answer
> 1 votes
The following may work, assuming you have an Illumina genome sequence library with decent coverage, preferably amplification-free:

Use kat to produce a kmer spectrum plot of the reads (kmer abundance vs. kmer frequency). When you look at, for example, the 21-mer spectrum of the reads, you would expect to see a peak in the spectrum where the kmer frequency corresponds with the estimated coverage of the genome. For instance, if your reads would be expected to provide 10x coverage of the original diploid genome, then you would expect a peak in the kmer spectrum at a kmer frequency around 10.

If the genome had been doubled, and each of the two copies had undergone a number of mutations, then you would also expect a second peak to be found with half the coverage of the first (at a kmer frequency of 5 in the example). The relationship between the heights of the two peaks would reflect the degree to which diploidisation had occurred post-doubling.

Also, I seem to remember something like a pachytene chromosome spread being used with an optical microscope.

---
Tags: cnv
---


# Question
Title: There's a faster way to download 1000 Genomes phase 3?
I'm trying to download the 1000 Genomes (phase 3) through FTP, but it's taking ages. I saw the documentation about using Aspera (fasp protocol), but the instructions don't work: `ERR [ascp] SSH authentication failed`. Eg:

```
ascp -i /home/ibmuser/.aspera/connect/etc/asperaweb_id_dsa.putty -Tr -Q -l 100M \
     -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi ./

```

Via Aspera Desktop client I got: `SSH_MSG_DISCONNECT: 2 Too many authentication failures`

Thanks!

# Answer
> 0 votes
You are using the command for older versions. In the page you link to, they give two commands, one for the newer and one for the older versions of aspera. You are using the one for the older. The only difference is the keyfile, but what you want is:

```
ascp -i /home/ibmuser/.aspera/connect/etc/asperaweb_id_dsa.openssh -Tr -Q -l 100M \
 -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi ./

```

At least, I just tried this a second ago with the aspera bundle I downloaded today and that's the command that worked: `asperaweb_id_dsa.openssh` and not `asperaweb_id_dsa.putty`.

---
Tags: 1000genomes
---


# Question
Title: Error with BWA Mem input having multiple fastq files using cat and process substitution
Running BWA Mem on a number of paired end fastq files using process substitution on the inputs results in this error:

```
[mem_sam_pe] paired reads have different names: "HWI-D00360:5:H814YADXX:1:2111:20284:33068", "HWI-D00360:5:H814YADXX:1:1213:17580:34847"

```

The command:

```
bwa mem -M ref.fa "<(cat U0a_CGATGT_L001_R1_004.fastq U0a_CGATGT_L001_R1_005.fastq U0a_CGATGT_L002_R1_001.fastq U0a_CGATGT_L002_R1_002.fastq U0a_CGATGT_L002_R1_003.fastq U0a_CGATGT_L001_R1_001.fastq U0a_CGATGT_L001_R1_002.fastq U0a_CGATGT_L002_R1_004.fastq U0a_CGATGT_L001_R1_003.fastq U0a_CGATGT_L002_R1_005.fastq)" "<(cat U0a_CGATGT_L001_R2_003.fastq U0a_CGATGT_L002_R2_005.fastq U0a_CGATGT_L001_R2_004.fastq U0a_CGATGT_L001_R2_005.fastq U0a_CGATGT_L002_R2_001.fastq U0a_CGATGT_L002_R2_002.fastq U0a_CGATGT_L001_R2_001.fastq U0a_CGATGT_L002_R2_003.fastq U0a_CGATGT_L001_R2_002.fastq U0a_CGATGT_L002_R2_004.fastq)" > out.sam

```

However the read (there were several) that threw this error does seem to be properly paired:

```
grep 'HWI-D00360:5:H814YADXX:1:2111:20284:33068' *.fastq
U0a_CGATGT_L001_R1_004.fastq:@HWI-D00360:5:H814YADXX:1:2111:20284:33068 1:N:0:CGATGT
U0a_CGATGT_L001_R2_004.fastq:@HWI-D00360:5:H814YADXX:1:2111:20284:33068 2:N:0:CGATGT

```

I can see that the concatenation of the input files is happening in a different order, so this is the next thing I'm going to try to fix, but has anyone else run into this?

(I probably just answered my own question but I'll post it anyway, that may or may not be the issue)

# Answer
> 3 votes
You are correct that the concatenation of input read files will confuse BWA. BWA depends on the input reads being in the same order in R1 and R2 files for paired-end alignment.

From bwa docs: `If mates.fq is present, this command assumes the i-th read in reads.fq and the i-th read in mates.fq constitute a read pair.` BWA will not search through FASTQ files trying to find matching mate IDs.

You will notice that your initial error message is reporting two different read IDs that are improperly paired. Presumably the mispaired read corresponds to the first read in the out-of-order file.

# Answer
> 3 votes
AS both the OP and @MaximilianPress suggested, **the paired-end reads should be provided as 2 separate arguments to bwa**. This is true whether the reads are in files or streamed. For the sake of completeness, I provide the corrected command here.

`bwa mem -M ref_file.fa <(cat A.read1.fq B.read1.fq ...) <(cat A.read2.fq B.read2.fq ...)`

**The order of the files is important**: if the order is A, B for read1, then it should be A, B and *not* B, A for read2.

# EDIT:

*(Thanks to @JohnMarshall for pointing this out.)*

**Avoid double quotes around process substitution. Quote the file names (if needed).** The quotes around process substitution, as apparently intended, work only on some systems/versions1, and seem like an undocumented feature to me.

**DON'T:**

`bwa mem -M ref_file.fa "<(cat A.read1.fq B.read1.fq ...)" "<(cat A.read2.fq B.read2.fq ...)"`

---

1 Oddly enough, process substitution inside double quotes works at least in the case I tested before posting: `bwa mem -M test1.fa "<(cat test_assembler_tiny_1.1.fastq)" "<(cat test_assembler_tiny_1.2.fastq)"` This contradicts the answer here: https://unix.stackexchange.com/q/486019 . I used bwa 0.7.17-r1188, with both GNU bash, version 3.2.57(1)-release (x86\_64-apple-darwin17) and zsh 5.7.1 (x86\_64-apple-darwin17.7.0), on macOS High Sierra 10.13.6. From OP's error msg, it seems it worked for @glyph too.

---
Tags: sequence-alignment, read-mapping, fastq, bwa, bash
---


title: HVG and variation coefficient

text:
# Question
Title: HVG and variation coefficient
I'm looking for highly variable genes from a data set. I have calculated the coefficient of variation of the expressed genes of my dataset and inserted them in a histogram graph. Now I should select only the highly variable genes. Should I select genes with a coefficient of variation \>1, or how? thank you for your help,

regards

# Answer
> 1 votes
Heuristically, I agree with @Jonathan Moore's answer. However, you could try to take a more disciplined approach and model the distribution as a mixture of two distributions, and infer the parameters of the two distributions from that, which will allow you to make somewhat more informed decisions. A basic intro to mixture modeling in R is here.

However, a worry of mine when I see data such as this is that you have some confounder that leads to that spike at CV ~ 1. I would actually plot the mean against the SD of the genes in question, and make sure that it is not due to observations of e.g. one transcript across samples when in most samples there are zero transcripts (see last section here). This would almost certainly be noise, and then what you would be enriching for is genes with very low expression, leading to a very small denominator for the CV calculation. If you are pseudocounting the expression (for example adding 1 to all values to avoid dividing by zero), I would try a few different values for the pseudocount to ensure that the value of the pseudocount is not significantly affecting your CV distribution that you show here. That would be a major warning sign of artifacts.

Therefore I would also suggest setting some sort of cutoff for minimal expression of genes in that distribution. Maybe an average of 1 transcript per observation (sample \[bulk RNA-seq\] or cell \[scRNA-seq\]) might be a reasonable cutoff for genes to consider.

I would do this data cleaning and filtering before making any decisions about this CV distribution.

# Answer
> 1 votes
There seems to be a dividing line around CoV=0.9. Perhaps this is a good place to start for a cutoff.

---
Tags: r
---


# Question
Title: How can I find all the columns available in UniProt.ws package for R?
I am trying to find the subcellular localization of my 10,000 proteins using UniProt.ws package for R. However, I am unable to find all the columns available for query.

I used another package named Bioconductor which had easy way to lookup all the columns available using COLUMNS(org.Hs.eg.db)

I tried reading the documentation and couldn't find anything which can output subcellular localization. I am new to this field, any help is appreciated.

# Answer
From the `UniProt.ws` manual:

> `columns` shows which kinds of data can be returned for the `UniProt.ws` object.

So, you will need to go over the output of `columns(your_uniprotws_object)` if it includes the data type you are looking for.

And here is a related question.
> 1 votes

# Answer
Do you have a small sample of the proteins you are looking for and what you have tried so far in R using Uniprot.ws ? (including output of it), it could be helpful to identify your problem and how to solve it.

Personally, I would rather consider Bioconductor as a platform for providing R packages dedicated to bioinformatics than a package to analyze data (but it's my opinion, maybe I'm wrong).
> 1 votes

---
Tags: r, bioconductor, uniprot, rstudio
---


title: Asking for Docking programme

text:
# Question
Title: Asking for Docking programme
I'm tring to do my M.Sc. research and I Have to do Docking to save some money as I cannot try all the compound that I'm working on.

what is the best Protein-Ligand Docking programme?

**And anyone has any tutorial for such programmes.**

# Answer
> 4 votes
# It's case specific

Docking small compounds well is problematic and you need to consider whether:

* if a single conformer is okay
* if rigid sidechains are okay
* if implicit water is okay
* if a non-polarisable model is okay
* if a non-covalent interaction is okay
* if you are okay with a score that is a rough ΔΔGibbs approximation of the interaction (see this SO answer) or do you want something more like a K\_d

Do note that for any experiment you do to have any validity you need to have controls. Do you know what can and cannot bind already? It is imperative that you score these too.

# Free

I am assuming that you want to use a free program. Virtual compound screens (VCS) are a big business in pharma. Therefore, a lot of really good software is not free, such as ICM-Dock, which outscores AutoDock discuss below in most tests.

## AutoDock: simple cases

AutoDock Vina is easy to use and has a PyMOL plugin (I believe), but the peptide, in both its sidechains and its backbone, is rigid, which is not good. The ligand rotamers are generated for AutoDock at runtime, but there is a hardcoded low limit, it technically isn't a rigid ligand, but nearly. It is the best-scoring free rigid peptide docking program.

There are lots of tutorials on using Autodock, such as the official one or on YouTube. As shown by this Bioinformatics Stack Exchange question Autodock 4 (different score fxn) has several steps, each important.

## Rosetta: very complicated cases

I prefer the harder to use Rosetta Ligand\_dock which uses conformers (if the ligand is parameterised correctly, say they were generate via Open Babel conformers) and allows sidechains to repack.

## Water

Both tools mentioned use implicit waters, which works badly for highly hydrophobic ligands. That is the waters are simulated as a homogeneous field that affects the ligand and protein. But some solutions are possible, in the case of AutoDock WaterDock variant exists —although this is not polarisable water. For more advanced cases there are MD methods such as dynamic undocking, which are completely non-trivial. All these methods mentioned use classical force-field calculations (Amber, Charmm, Talaris etc. are different models), but in MD noise (think of it as environmental heat) is present in the calculations, so the ligand comes unbound (and is actually pulled out). You can also get better results by refining your bound ligand with QM-MD (*e.g.* gauss of gaussian).

## Choices...

This is a lot of options. So a first question is are your ligands man-made drugs/fragments (i.e. obey Lipinski's rule of 5) or metabolites? Generally the former are rigid and very hydrophobic (greasy) (and require better force fields), while the latter are flexible and hydrophilic.

## Minor considerations

Also to consider is your entranceway to your protein as the Michaelis constant is dictated in some enzymes by the substrate going in, not by its binding (e.g. P450s).

If the ligand have a shared backbone that is fixed there are some tricks you can do to cut your compute time and increase the reliability.

If the ligands are bound not all methods are possible. Rosetta ligand\_dock can but it is a bit fiddly to parameterize.

## RDKit skills

One thing also to consider is your programming skills. Can you use RDKit? In most cases downloading the sdf off PubChem and running OBabel or similar to generate conformers is not possible and you have be familiar with a python package called RDKit, which is really powerful but non-trivial to use.

Additionally, docking scores can be strongly augmented with ligand-only values (such as logP, TPSA, molar refractivity, molecular weight, QED, Bertz topology etc. etc) that are generated from RDKit.

NB. The scoring function in Autodock Vina (not 4 though, as that is Amber FF) partially takes into account.

## Ligand\_dock details

Ligand\_dock is an older application of the Rosetta suite and the current philosophy is to use a specific Rosetta script, however, it works fine and is actually quicker to implement. Here are some links:

# Answer
> 2 votes
`Autodock-Vina` is a free molecular docking solution. You need to generate different input conformations for each ligand though, e.g. with RDKit.

Here is a link to tutorial

---
Tags: docking
---


title: Post docking evaluation

text:
# Question
Title: Post docking evaluation
What to do after protein-ligand docking? Is there any necessary or enhancement step recommended after docking to validate the confirmation of highest score is the best? For your information, I am using CCDC GOLD for docking and GOLD score as my scoring function.

# Answer
You can run a Molecular Dynamics simulation to see if the pose is stable and whether the protein changes/adapts it’s conformation.

You can calculate the ligand efficiency, the docking score divided by the number if heavy atoms, and rank the compounds accordingly.

You can also compare the scores of known ligands to the experimental values and see if the docking the scores are predicting known ligands or correlate with the experimental ones.

At some point you have to make a selection and purchase or synthesize the compounds. How many compounds can you test experimentally? How many can you purchase? What is the purpose of these ligands going to be? Drugs? Did you screen the correct set of ligands e.g. a library of lead compounds that are ready to purchase for cheap, that would be a good start.
> 2 votes

---
Tags: proteins, protein-structure, docking, protein-protein-interaction
---


# Question
Title: how to calculate Pairwise Alignment Scores for blosum62
I have a bioinformatics exam coming up. I can understand the difference between the global algorithm and the local algorithm, but I have a problem with gap opening penalty and gap extension penalty. What is the difference between gap opening penalty and gap extension?

Can anybody help me in working out how to use the gap opening penalty and gap extension penalty for blosum62?

There are a lot of questions like the one below:

*Compute the score of the following sequence alignment given the BLOSUM62 matrix below and gap opening penalty (GOP) = 12, and gap extension penalty (GEP) = 2.*

```
SDRVIKAAIFDPIQPDF---G-----PVYFGLGHVH
RDLVERLFILDMI-PGLIKAGDSFPIPVALMINHIF

   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V
A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0
R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3
N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3
D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3
C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1
Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2
E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2
G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3
H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3
I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3
L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1
K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2
M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1
F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1
P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2
S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2
T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0
W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3
Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1
V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4

```

# Answer
Your alignment is

`SDRVIKAAIFDPIQPDF---G-----PVYFGLGHVH
RDLVERLFILDMI-PGLIKAGDSFPIPVALMINHIF`

For each position in the alignment you calculate the score for that alignment. For position 1 we'd look up S vs R in the matrix and find a score of -1. You continue doing this until you hit the first -, which is not in the matrix. The first - is a gapopening, each consequent - in a series of -'s counts as a gap extension, instead of an opening.

`---G--
IKAGDS` gives `-12 -2 -2 6 -12 -2`
> 1 votes

---
Tags: sequence-alignment, homework
---


# Question
Title: How I can change the name of multiple files at once in R or terminal?
I have 200 .vcf files in a folder with long names like

```
LP6005409-DNA_E03_vs_LP6005408-DNA_E03.snp.pass.vcf

```

How can I change the name of each file to, for example `A1.vcf ....A200.vcf`?

Renaming these files manually takes ages, help please.

```
> file.rename(list.files(pattern = "vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf"))
Error in file.rename(list.files(pattern = "vcf"), paste0("A", 1:length(list.files(pattern = ".vcf")),  :
  'from' and 'to' are of different lengths

[fi1d18@cyan01 new]$ for i in *.vcf; do echo mv $i A${k}.vcf ; k=`echo "1+${k}" | bc`; done
mv LP2000325-DNA_A01_vs_LP2000317-DNA_A01.snp.pass.vcf A.vcf
(standard_in) 2: syntax error
mv LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf A.vcf
(standard_in) 2: syntax error
[fi1d18@cyan01 new]$ for i in *.vcf; do mv $i A${k}.vcf ; k=`echo "1+${k}" | bc`; done
(standard_in) 2: syntax error
mv: try to overwrite `A.vcf', overriding mode 0500 (r-x------)?

```

# Answer
> 5 votes
just to complete your question, you can do it in R in one line.

After setting the directory as your working directory, just type the following:

```
file.rename(list.files(pattern = "vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf"))

```

Maybe you can try `file.rename(list.files(pattern = ".vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf"))`. Also, if there is still an error, can you post the output of `list.files(pattern = ".vcf")` just to see the structure of your files in folder (maybe one file has a weird name not detected by the pattern).

# Answer
> 9 votes
Use `perl-rename`. This is usually called `rename` on Debian-based systems like Ubuntu or Mint, and `perl-rename` or `prename` on others. Assuming you have it as `rename`, you can simply do:

```
rename -n 's/.*\.vcf/"A" . ++$c . ".vcf"/e' *snp.pass.vcf

```

The `-n` causes `rename` to only print what it would do, without doing anything. So once you make sure it does what you want, run the command again without the `-n`.

# Answer
> 3 votes
If your filenames do not contain blank spaces, you can do it with a for-loop in `bash`:

Make a test run with `echo` so that the command is only printed but not yet executed:

```
for i in *.vcf; do echo mv $i A${k}.vcf ; let k++ ; done

```

And if the result looks good remove the `echo` statement to execute the `mv` commands:

```
for i in *.vcf; do mv $i A${k}.vcf ; let k++ ; done

```

# Answer
> 3 votes
**EDIT:**

*(Based on the comments by @terdon, with minor changes)*

Try this Perl one-liner. It is recommended for general use, even if your file names contain "unexpected" characters such as newlines.

```
perl -e 'for ( 0..$#ARGV ) { rename $ARGV[$_] => "A@{[ $_ + 1 ]}.vcf"; } ' *pass.vcf

```

Here, `-e` tells `perl` to look for code on the command line instead of a file with the script. `@ARGV` is the array of the command line arguments, here, files with old names.

`for ( 0..$#ARGV ) { ... }` loop goes through the indexes of elements of this array, assigning each implicitly to `$_`. The first argument to `rename` is the old file name: `$ARGV[$_]`. Instead of the comma separator of the method arguments I use `=>` synonym for clarity in this method call.

The new file name is the second argument. `@{[ ... ]}` inside the double quotes is the Perl idiom to evaluate the expression inside, which is `$_ + 1`, to make the file numbering start with 1. If you don't care about this, it is easier to make it start with 0: just use `"A$_.vcf"`, or, better, `"A${_}.vcf"`. The last form needs to be used for cases like this: `"A${_}pass.vcf"` to prevent Perl from treating `$_pass` as a variable, where `$_` is intended.

**DEPRECATED:**

Use the one-liner below when you are sure that the old file names do not have "unexpected" characters such as newlines. It is not for general use, see: Why you shouldn't parse the output of ls.

It will work with the example shown by the OP. It also works with blanks and tabs in the file names, but it **breaks when file names have newlines**.

The (deprecated) Perl one-liner:

```
ls *.pass.vcf | perl -lne 'rename $_, "A$..vcf";'

```

Here, `ls |` lists all files with old names, one per line and passes the stdout into stdin of the next command, which is the Perl one line script. `-l` removes the new lines from the input. `-n` reads the input one line at a time. These are all commonly used command line Perl options.

In the one liner script, `$_` is the line read from stdin, and `$.` is the line number.

---
Tags: r, perl, linux, shell, awk
---


# Question
Title: What is a samtools mpileup reference skip?
The samtools documentation for mpileup states:

> At this column, a dot stands for a match to the reference base on the forward strand, a comma for a match on the reverse strand, a '\>' or '\<' for a reference skip ... Similarly, a pattern \`-\[0-9\]+\[ACGTNacgtn\]+' represents a deletion from the reference.

What is a reference skip?

# Answer
That documentation has been expanded recently (see PR #1055) and now describes the characters seen in this column in more detail:

```
Forward   Reverse                    Meaning
---------------------------------------------------------------
 . dot    , comma   Base matches the reference base
 ACGTN     acgtn    Base is a mismatch to the reference base
   >         <      Reference skip (due to CIGAR "N")
   *         *      Deletion of the reference base (CIGAR "D")

```

(And so on — see the pull request or *doc/samtools-mpileup.1* on the **develop** branch for the complete pileup format description. The documentation on that web page will be updated when this is released in a new samtools release.)

As noted in https://bioinformatics.stackexchange.com/a/421/134, SAM has always had separate `D` (deletion from the reference) and `N` (skipped region from the reference) operators. These are treated similarly computationally but emitted intentionally by aligners to distinguish introns and the like from other deletions.
> 5 votes

# Answer
Skipped references are similar to deletions, but the different symbols used are indicating that it's an *expected* deletion. The most common use for a skipped reference is when excluding intronic sequence from cDNA reads that are being matched to a genome reference.

Note that the skipped reference bases are added by the mapping program, which generates the SAM/BAM/CRAM format and decides on whether to call a deletion or a skipped base.
> 1 votes

---
Tags: samtools, mpileup
---


# Question
Title: How to loop multiple function in shell script?
I need to extract sequences one after another consecutively from a large fasta files (multiple fasta files) and each extracted files to be saved in new fasta file (I mean the first sequence extracted from multiple fasta file to be saved as gene1 and second sequence extracted from multiple fasta file to be saved in another file and so on). For that, I tried a script like this,

```
for file in *.fasta do cat "$file" \
  | awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} \
          {printf("%s",$0);} \
     END {printf("\n");}' \
  | awk '{if ( NR==1 ) print $0;}' \
  | awk '{print $1"\n"$2}' | uniq > "new_${file}"

```

But, its not serving my purpose. Firstly, I could not loop my sequence extraction procedure from multiple fasta file and save as gene1 and gene2 and so on.

```
For example,

```

I have fasta files as follows,

```
org1.fasta
>locus1
ATGCGTAGAG
>LG12
TATGCATGAT
>bpl
TAGATGGAT

org2.fasta
>ysr
ATGTAGCGA
>inc
TGATCTATG
>pym
ATGCATAGA

org3.fasta
>siv
TAGTAGTAT
>tdp
ATGAGTCGA
>tem
ATGCTHATG

```

I need to extract first positioned sequences from all the `org1.fasta org2.fasta and org3.fasta` files and save in a new file `gene1.fasta` and second positioned sequences as `gene2.fasta`and third positioned sequences to be saved as `gene3.fasta`.

```
Expected output,

gene1.fasta
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT

gene2.fasta
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA

gene3.fasta
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

```

Please help me to do the same. Thanks in advance.

# Answer
You don't need a loop. You can do the whole thing with a simple awk one-liner:

```
awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta

```

I ran this on your files and got:

```
$ for file in gene*; do echo "=== File: $file ==="; cat $file;  done
=== File: gene1.fasta ===
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT
=== File: gene2.fasta ===
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA
=== File: gene3.fasta ===
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

```

### Explanation

* `-vRS='>'` : set the record separator to `>`. This will make `awk` consider every fasta sequence as a single "line".
* `FNR>1{foo}` : `FNR` is the line number of the current file. Because `>` is the first character of the entry, the first "line" will always be empty, so it needs t obe skipped.
* `printf ">%s",$0 > "gene"FNR-1".fasta"` : this does the work. It prints a `>` followed by the current "line" (`$0`, this is the entire sequence + the header) and redirects it into a file called `gene`, then the current line number -1 (-1 is needed because of the issue mentioned above) and a `.fasta` extension.

=====

Note that `org*fasta` will be expanded with the default collation order. This means that `org10.fasta` will be processed before `org1.fasta`:

```
$ printf "org1.fasta\norg10.fasta\n" | sort
org10.fasta
org1.fasta

```

To avoid this and have the files in numerical order, you can use brace expansion instead of the `*`. So, if you know your files are all `orgN.fasta` where `N` can take values from 1 to 20, do this:

```
awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org{1..20}.fasta

```
> 3 votes

# Answer
The title of your question indicates a `shell` solution, however you have added `R` and `Python` as tags. Here is an `R` solution making use of `Biostrings` pacakge (using two input files for simplicity):

```
library(Biostrings)

# create Biostring objects from each multifasta file
org1 <- readDNAStringSet("org1.fasta")
org2 <- readDNAStringSet("org2.fasta")

# loop over the Biostring objects
for(i in 1:length(org1)){

  output_file <- paste0("gene_", i, ".fasta")

  # create a Biostring object for each index
  genes_at_index_i <- c(org1[i], org2[i])
  assign(gene_name, genes_at_index_i)

  # save Biostring object as a multifasta file
  writeXStringSet(genes_at_index_i, output_file)
}

```

Output, gene 1:

```
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA

```

Output, gene 2:

```
>LG12
TATGCATGAT
>inc
TGATCTATG

```

Warning: The code above won't work if you don't have the same number of fasta entries in each multifasta file.
> 1 votes

---
Tags: r, python, bash, perl, shell
---


# Question
Title: WGCNA module preservation error code with plotting summary statistics
I am trying module preservation with WGCNA using my own data. I am plotting the summary statistics including Zsummary, medianRank, observed IGP and IGP permutation p-value for two comparison data sets.

When plotting the scatterplots for IGP & p-value vs. module size I get following error code. I have re-run code and done practice data with no problems but unsure how to correct for this error? Does anyone have any experience with this and know how to resolve? Thank you kindly

```
 plotData[[2]][, 3] = cr[[ref]][[test]]$Acutal.IGP
  plotData[[2]],[, 4] = -log(cr[[ref]][[test]]$p.value + 1e -4)

errror in plotData[[2]][, 3] = cr[[ref]][[test]]$Acutal.IGP
number of items to replace is not a multiple of replacement length

```

# Answer
> 1 votes
My first guess is that there is a typo in your code, and Acutal.IGP should be Actual.IGP. I checked my tutorials and the word Acutal is in none of them, so I would assume you mistyped the variable name.

---
Tags: r, wgcna
---


title: AutoDock 4 -soft docking

text:
# Question
Title: AutoDock 4 -soft docking
How could I turn on/off the smooth potential in AutoDock4? I need to do two docking. The first one is soft docking with smooth potential and the second is re-docking without the smooth potential.

# Answer
> 1 votes
The `smooth` parameter is an optional keyword and command that changes the smoothing potential. If not provided, the default is 0.5 angstroms by default. Try passing 0 with this keyword in the AutoGrid Parameter File to turn off smoothing potential.

See the AutoDock 4.2.5 documentation for more details.

---
Tags: proteins, software-usage, docking
---


# Question
Title: How to merge two files on the first 3 columns?
I have two files that look like this:

```
$ head ret
chr   pos     gene_id      pval_ret
chr1 10177 ENSG00000227232 0.381708
chr1 10352 ENSG00000227232 0.959523
chr1 11008 ENSG00000227232 0.218132
chr1 11012 ENSG00000227232 0.218132
chr1 13110 ENSG00000227232 0.998262
chr1 13116 ENSG00000227232 0.438572
chr1 13118 ENSG00000227232 0.438572
chr1 13273 ENSG00000227232 0.694523
chr1 14464 ENSG00000227232 0.825113

```

```
$ head lcls
chr  pos     gene_id       pval_nominal
chr1 13550 ENSG00000227232 0.375614
chr1 14671 ENSG00000227232 0.474708
chr1 14677 ENSG00000227232 0.699887
chr1 16841 ENSG00000227232 0.127895
chr1 16856 ENSG00000227232 0.627822
chr1 17005 ENSG00000227232 0.802803
chr1 17147 ENSG00000227232 0.443987
chr1 17407 ENSG00000227232 0.397521
chr1 17408 ENSG00000227232 0.615126

```

How can I merge these two files, `ret` and `lcls`, on the first 3 columns: chr, pos, gene\_id?

I know how to merge them on the first two columns:

```
awk 'NR==FNR{a[$1,$2]=$3;next} ($1,$2) in a{print $0, a[$1,$2]}' ret lcls > overlap

```

But how can I do it on the first 3 columns? In the output I want these 4 columns: `chr`, `pos,gene_id`, `pval_ret`, `pval_nominal`.

# Answer
> 2 votes
You can use the same approach as for two fields, just add the third:

```
awk 'NR==FNR{a[$1,$2,$3]=$4;next} ($1,$2,$3) in a{print $0, a[$1,$2,$3]}' ret lcls

```

# Answer
> 0 votes
With the 'paste' command you can merge two files. Example :

```
paste ret lcls | awk '{print $1" "$2" "$3" "$6" "$4" "$8}'

```

Do you want 4 or 5 columns ? A sample of the result you want could help understand your needs.

---
Tags: text-processing, awk
---


# Question
Title: different gene expression in the monocle
As it is shown in the picture that we can get the pseudotime plot by using Monocle. I was wondering whether we can use some function like the Seurat FindMarker to find out the differential expression genes between the Double-KO and WT in the bottom right corner?

# Answer
Seurat FindMarker Doc Yes based upon the R documentation here it seems possible.

This Differential expression article seems to show this is how.

```
library(Seurat)
pbmc <- readRDS(file = "../data/pbmc3k_final.rds")
# list options for groups to perform differential expression on
levels(pbmc)
# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
monocyte.de.markers <- FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono")
# view results
head(monocyte.de.markers)
# Find differentially expressed features between CD14+ Monocytes and all other cells, only
# search for positive markers
monocyte.de.markers <- FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = NULL, only.pos = TRUE)
# view results
head(monocyte.de.markers)
# Pre-filter features that are detected at <50% frequency in either CD14+ Monocytes or FCGR3A+
# Monocytes
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", min.pct = 0.5))
#...
# Test for DE features using the MAST package
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", test.use = "MAST"))
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))

```

Maybe something like this would work for you. Monocle export

```
  #replace the monocle_cds with your monocle
  seurat <-exportCDS(monocle_cds, export_to = c("Seurat", "Scater"))
  #This bellow will list the options for ident.1 and ident.2
  levels(seurat)
  # insert name from levels(seurat) command in parentheses
  head(FindMarkers(seurat, ident.1 = "Double-KO+ Mono", ident.2 = "WT+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))

```

Hopefully you are able to separate the two groups Double-KO and WT into ident.1 and ident.2 so you can use the function on them. Good luck. As long as the exportCDS function works on your data type you input for cds.
> 1 votes

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


# Question
Title: GenomicRanges: get nearest neighbor distances for random genes using a for loop
I want to compute the nearest neighbour distances using a for loop. I can do this with a random subset by doing this:

```
df <-  read.table("tmp.txt", header=F)
colnames(df) <- c("gene", "seqnames", "start", "stop", "LOCID")
gr <- makeGRangesFromDataFrame(df, keep.extra.columns=TRUE)

# genes of interest in a list
gene.list <- read.table("careplastic_diffs.txt")

# filters the GenomicRanges object with genes of interest
sub.gr <- gr[(elementMetadata(gr)[, "gene"] %in% gene.list$V1)]

# compute distances
dists <- distanceToNearest(x = sub.gr)

mcols(dists)[,1] # contains a list of distances

```

I now want to essentially do this 5 times, using the same logic as above:

```
## create a rand dataframe, each column is a random list of 50 genes
rand.df <- data.frame(replicate(5, sample(gene.list$V1, 50)))

# loop through columns and gets distances and prints to a file
for(i in 1:ncol(rand.df)){
  rand.gr[i] <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists[i] <- distanceToNearest(x = rand.gr[i])
  output <- print(mcols(dists[i][,1]))
}

```

But I can't seem to make this work.

> Error: subscript contains out-of-bounds indices In addition: Warning message: In recycleSingleBracketReplacementValue(value, x, nsbs) :
> number of values supplied is not a sub-multiple of the number of values to be replaced

It's probably something really silly but any help would be appreciated ...

# Answer
Is it because rand.gr has not been initialised, so when you ask to set rand.gr\[i\] it can't because there isn't one yet?

If you set i to a value (e.g. 1), and run the lines inside your loop, which of the lines produces the error?

If you don't need to retain all the random genomic ranges at the end, you could just leave out some of the \[i\]s, something like this:

```
for(i in 1:ncol(rand.df)){
  rand.gr <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists <- distanceToNearest(x = rand.gr)
  output <- print(mcols(dists[,1]))
}

```
> 1 votes

# Answer
You could do this much more easily and much faster on the command line, starting with an unsorted BED file of genes called `genes.unsorted.bed`:

```
$ sort-bed genes.unsorted.bed | shuf - | head -5 - | sort-bed - | closest-features --closest --dist --no-overlaps - > answer.bed

```

Replace the `5` in `head -5` with however many random genes you want. Remove `--no-overlaps` if you want to allow consideration of overlapping genes (which returns a distance of zero).

Output columns are delimited with `|`. This can be replaced with another delimiter, but the output can be brought into R easily with `read.table` or similar.

Ref: https://bedops.readthedocs.io/en/latest/content/reference/set-operations/closest-features.html
> 0 votes

---
Tags: r, bioconductor, genomicranges, genomefeatures
---


# Question
Title: In silico PCR: any models that can quantitatively predict the initial content from the end amplified result?
Let's say I've run the PCR and I now have a certain amount of the amplicon product. Is there any model that infers how much of it was there before the PCR started? I suspect it should be a probabilistic model (a distribution) parameters (shape) of which will strongly depend on the details of the PCR protocol and the characteristics of the amplicon...

If not, then may be there are forward models that predict the amount of the amplified result from the initial content?

# Answer
> 1 votes
The basic equation for PCR kinetics (Equation 1) states that the amount of amplicon after c cycles (Nc) is the starting concentration of the amplicon (N0) times the amplification efficiency (E) to the power c. The PCR efficiency in this equation is a number between 1 and 2 (2 indicates 100% efficiency). NCBI

In this article on research gate they said amplifx was the software you were looking for to, calculate "(i.e. expected efficiency in PCR reaction)". The amplifx even mentions "Testing primers in in silico PCR".

---
Tags: pcr
---


# Question
Title: How I move these files based on their names in R or terminal
> EDITED

I have list of my samples in first column and corresponding file name for my .vcf files in second column of a file name `clin_name.txt` like below. For example name of one vcf is `LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf`

```
OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

```

This is the link of my file https://www.dropbox.com/s/mvzf47feiqe7x69/clin\_name.txt.txt?dl=0

In the second file name `pre_post.list` I have information of pre or post treatment for each sample like this

```
pre-treated post-treated
OC/AH/120   OC/AH/065
OC/AH/126   OC/AH/069
OC/AH/136   OC/AH/074
OC/AH/160   OC/AH/084

```

This is link of this https://www.dropbox.com/s/ikrl1h8yoc0sngx/pre\_post.list.xlsx?dl=0

As you suggested I typed

```
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

```

Then I have

```
pre-treated
120
126
136
160
182

post-treated
065
069
074
084

```

in I typed `./test.sh` in which I wrote

```
for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/
    else echo post_list $item;  fi;  done;

```

That returns

```
cp: cannot stat `*OC/AH/069/S/T3*': No such file or directory
[fi1d18@cyan01 snp]$

```

Or I run

```
for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/pre/
    else cp *$out* /temp/hgig/fi1d18/TRG45/snp/post/;  fi;  done;

```

That returns

```
cp: cannot stat `OC/ST/026/S/T1*': No such file or directory
[fi1d18@cyan01 snp]$

```

I really don't know how to solve this :(

# Answer
> 1 votes
Assuming the names of all the files from one group are in file named list.

```
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

```

Now you have two number list

```
    for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_pre_files/
    else echo post_list $item;  fi;  done

```

for post

```
    for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_post_file/
    else echo Prelist $item;  fi;  done

```

or you could do this if one that isn't in post belongs in pre

```
    for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_post_file/
    else cp *$out* /dir/for_pre_file/;  fi;  done

```

---
Tags: r, bash, looper
---


title: Peak-calling using homer

text:
# Question
Title: Peak-calling using homer
I have a total of 78 ATAC seq samples from which Im trying to do a peak call .I tried this `batchParallel.pl findPeaks peaks` function but i couldn't find the output what it makes.

So I tried using loop but again the peak file is not written due to incorrect argument usage I guess.

These are my sample directories where i have all the tags

```
SRR2920564  SRR2920565

```

I used homer to make tag directories and now i using it for peakcalling.

here is my command.

```
for d in */ ; do findPeaks $d -style factor -o "${d}.txt" ; done

```

What i need is all the peakfile should be written as the directory name.

Such as if directory is "`SRR2920564`" then my peak file should be written in the same directory as `SRR2920564.txt`

I do run the loop it runs successfully but no output written. May be a mundane mistake am not able to figure it out.

Any suggestion or help would be highly appreciated .

# Answer
You just need to strip the `/` from `$d`, for which there are a number of options:

```
for d in */ ; do
    dname=`basename $d`
    findPeaks $d -style factor -o ${dname}.txt
done

```

or

```
for d in */ ; do
    dname=${d%%/}
    findPeaks $d -style factor -o ${dname}.txt
done

```

There are likely other ways one could go about this. Having said that, since you're not running these in parallel you could simply move the default output file (`peak.txt`) to a file name of your choice (using one of the strategies above).
> 1 votes

---
Tags: shell, atac-seq, peak-calling
---


# Question
Title: ATAC-seq peak annotation and downstream analysis
I have 7 different cell types, 4 normal cell types and 3 disease cell types. I have done until peak calling. Now I want to do peak annotation.

The sample I have is all biological replicates. Should I merge the peak file and then do peak annotation, or should I be doing individual peak annotation?

# Answer
The way we do this (and I think its quite common?), is to merge all the files of the same condition (taking the same number of reads from each), and then calling peaks on the merged sample - so you would have two peak sets - one from normal and one from disease. You then do the union of the two peak sets.

This makes annotation much easier.

For downstream quantitative analysis, you can then quantitate each individual sample against the merged peak set.
> 1 votes

---
Tags: atac-seq, homer
---


# Question
Title: Bedtools get fasta and ORF from a blastX run
Hello I made a `blastX research` (A query genome in nucleotide format translated into protein in 6 reading frames against a protein dabatase)

And here is a head of the result :

```
IDBA_scaffold_7517  493 YP_009316100    0.360   61  27  1   150 4   126 186 6.775E-03   37
IDBA_scaffold_29149 519 AIG51500    0.428   35  18  1   2792    2896    392 424 2.712E-03   40

```

I have 2 hits:

The first one is the query `IDBA_scaffold_7517` that have a hit with a protein sequence `YP_009316100` in the `(-)` strand because the start begins by `150` and ends with `4`.

The seconde one is the query `IDBA_scaffold_29149` that have a hit with a protein sequence `AIG51500` in the `(+)` strand because the start `begins 2794` and `ends by 2893`.

so with that I actually created a bed file in order to extract the fasta sequences from my genome ( I assume that since there is a hit with a protein in the Genomes the coordinates of this hit should be in a good ORF (without stop codon).

Here is the `bed_file` I made :

```
IDBA_scaffold_7517  4   150 Species 0   -
IDBA_scaffold_29149 2794    2893    Species 0   +

```

Then in order to extract the sequences from the coordinates in the Genome I used:

```
bedtools getfasta -fi Genome.fa -bed bed_file -s -fo test

```

with `-s` to force strandedness. If the feature occupies the antisense, strand, the sequence will be reverse complemented.

But when I do that the sequences extracted from the query are not in the good ORF. I guess it is because I should take into account the fact that the coordinates are from nuc translated to portein but how could I do that? Should I change the coordinates manually from the bed file according to the fact that it is a negative of positif strand ?

Thank you for your help.

# Answer
> 2 votes
Blast uses a 1-based coordinate system. Bed is 0-based half open. Your code needs to convert appropriately. See for example: https://www.biostars.org/p/317112/ .

---
Tags: fasta, bedtools, orf
---


# Question
Title: What are the standard ways to visualize protein-protein or gene-gene interactions?
I would like to visualize this interaction list. Is there an online/web based way to do it? Is there a way to analyze the data?

Assuming I have exported the list into lines such as:

`2 A2M 348 APOE`

# Answer
since I am not yet allowed to comment, I will have to pots this as an answer.

I guess Cytoscape (http://cytoscape.org/) would be a nice way to analyse such interaction lists/networks (still depends a bit on what you plan to do).

However, for this you will have to preprocess the downloadable excel sheet a bit.

1. save as tab-separated file (tsv)
2. run the below python3 script to format into an easily importable format for cytoscape. You can run the script as `python3 script.py <path-to-tsv-file> <delimeter for cytoscape>` . The new format is written to std out on the console, therefore maybe you need to write the output to a new file. A sample call could look like `python3 testbench.py human_interactome_2017-06-12.tsv.csv ";"`. Make sure to add the delimeter in quotes.
3. you can then import the new file into cytoscape (import network as table) where the first column is your source, the third column your target and the middle column your interaction type.
4. Within cytoscape you can then adapt your edge style to vary the line type in a discrete mapping depending on the value of the interaction.

Screenshot from cytoscape:

Code below:

```
import sys

HIIIidx = 4
HII5idx = 5
venkatesanidx = 6
yuidx = 7
Litidx = 8

print(sys.argv)

if len(sys.argv) <= 1:
    print("Usage: script.py file delimeter")

delim = ';'
inFile = sys.argv[1]

if len(sys.argv) > 2:
    delim = sys.argv[2]

def printInteraction(syma, symb, inter):
    print(syma + delim + inter + delim + symb)

with open( inFile , 'r') as file:

    file.readline()

    for line in file:

        aline = line.strip('\n').split('\t')

        if aline[HIIIidx] == '1':
            printInteraction(aline[1], aline[3], 'HIII')

        if aline[HII5idx] == '1':
            printInteraction(aline[1], aline[3], 'HII5')

        if aline[venkatesanidx] == '1':
            printInteraction(aline[1], aline[3], 'VEN')

        if aline[yuidx] == '1':
            printInteraction(aline[1], aline[3], 'YU')

        if aline[Litidx] == '1':
            printInteraction(aline[1], aline[3], 'LIT')

```
> 2 votes

# Answer
Your list has two identifiers for the same node per line. In order to use it, you will need to change that. If you want to use the gene name (2nd and 4th fields, in your example), just run:

```
awk 'print $2,$4' netw.txt > netw.gr

```

If you want to use the Entrez geneIDs instead, run:

```
awk 'print $1,$3' netw.txt > netw.gr

```

Then, as others already mentioned, install `Cytoscape`, launch it and import your gene list:

1. In the welcome screen, select "From Network File...":
2. Then, select your network file (`netw.gr` in the example above) and choose "Advanced Options":
3. Set the delimiter to SPACE and uncheck "Use first line as column names":
4. Finally, click on the header of each column and set the first to "source node" and the second to "target node":

Now, click "OK" and you will have successfully imported your network into Cytoscape.
> 2 votes

# Answer
Usually interactions are represented in a graph, with edges representing the interactions (colors, thickness or values on the edges represent different properties of those edges).

As far as I know, there aren't on line programs to represent this kind of information. For an off line representation, you can use the igraph package in R, or Cytoscape program.
> 1 votes

# Answer
Several have mentioned the excellent Cytoscape desktop package for visualising interactions as a graph. If you want to build your own web-based visualisations (rather than using a pre-built website to visualise them) then the Cytoscape javascript library might fulfil your needs.
> 1 votes

---
Tags: gene, proteins, interactions
---


title: Binding columns with matched elements

text:
# Question
Title: Binding columns with matched elements
I have a matrix in which sample IDs are in rows and genes are in column like below

```
> head(data)
                                 ETV7    PLAUR  SLC11A1
TCGA.JY.A6FB.01A.11R.A336.31 1.791167 6.241607 3.709848
TCGA.JY.A938.01A.11R.A37I.31 5.389175 7.048782 3.124435
                                BIRC3 TNFRSF1B
TCGA.JY.A6FB.01A.11R.A336.31 6.726991 5.758026
TCGA.JY.A938.01A.11R.A37I.31 8.097783 6.621900
>

```

In another matrix I have sample IDs, months to death and status while the order of sample IDs is not the same with my previous matrix like below

```
Patient_ids OS_time  OS_event
TCGA-IG-A4QS    3.933333333 1
TCGA-JY-A6FB    61.23333333 0
TCGA-JY-A938    35.33333333 0
TCGA-JY-A939    22  0
TCGA-JY-A93C    23.5    0
TCGA-JY-A93E    25.56666667 0
TCGA-L5-A4OJ    21.3    0
TCGA-L5-A4OS    59.4    0
TCGA-L5-A8NM    7.866666667 1
TCGA-R6-A6KZ    5.133333333 1
TCGA-R6-A6XG    38.93333333 0

```

How I can bind OS\_time and OS\_event columns from second matrix to the first matrix keeping the right sample order?

---

> EDITED

---

```
> # Reproducible example:
> patient_id = c("TCGA.JY.A6FB.01A.11R.A336.31","TCGA.JY.A938.01A.11R.A37I.31","TCGA.L5.A4OS.01A.11R.A212.31","TCGA.AA.A400.01A.11R.A212.31")
> ETV7 = c(1.79,5.38,5,6)
> PLAUR = c(6.24,7.05,5,6)
> SLC11A1 = c(3.7,3.12,5,6)
> df = data.frame(ETV7,PLAUR, SLC11A1)
> rownames(df) = patient_id
> # OS data
> Patient_ids = c("TCGA-IG-A4QS","TCGA-JY-A6FB","TCGA-JY-A938","TCGA-JY-A939","TCGA-JY-A93C","TCGA-JY-A93E","TCGA-L5-A4OJ","TCGA-L5-A4OS")
> OS_time = c(3.93,61.23,35.33,22,23.5,25.56,21.3,59.4)
> os_data = data.frame(Patient_ids,OS_time)
> library(rebus.base)
> df$Patient_id = gsub(DOT,"-",rownames(df))
> position = sapply(os_data$Patient_ids,function(x){grep(x,df$Patient_id)})
> position = lapply(position,function(x){c(x,rep(NA,1-length(x)))})
> os_data$Position_matched = unlist(position)
> m = data.frame(matrix(ncol = length(df),nrow = nrow(os_data)))
> colnames(m) = colnames(df)
> for(i in 1:length(position))
+ {
+   if(is.na(position[i])){}
+   else{
+     m[i,] = df[position[i],]
+   }
+ }
Error in xj[i] : invalid subscript type 'list'
>

```

# Answer
Here a little way to do it. I don't think it is a perfect solution but it should work. The trick here is to compare `patients_id` that are wrote with a `.` or with a `-`. So, here is how I do (I modified your example in order to reflect several configurations):

```
# Reproducible example:
patient_id = c("TCGA.JY.A6FB.01A.11R.A336.31","TCGA.JY.A938.01A.11R.A37I.31","TCGA.L5.A4OS.01A.11R.A212.31","TCGA.AA.A400.01A.11R.A212.31")
ETV7 = c(1.79,5.38,5,6)
PLAUR = c(6.24,7.05,5,6)
SLC11A1 = c(3.7,3.12,5,6)

df = data.frame(ETV7,PLAUR, SLC11A1)
rownames(df) = patient_id

# OS data
Patient_ids = c("TCGA-IG-A4QS","TCGA-JY-A6FB","TCGA-JY-A938","TCGA-JY-A939","TCGA-JY-A93C","TCGA-JY-A93E","TCGA-L5-A4OJ","TCGA-L5-A4OS")
OS_time = c(3.93,61.23,35.33,22,23.5,25.56,21.3,59.4)
os_data = data.frame(Patient_ids,OS_time)

```

Now, I create a new column in `df` that will contains `patient_id` wrote with a `-`:

```
library(rebus.base)
df$Patient_id = gsub(DOT,"-",rownames(df))

```

Now, I'm looking for common pattern between `df$Patient_id` and `os_data$Patient_ids`, find their matching row position and extract these lines:

```
position = sapply(os_data$Patient_ids,function(x){grep(x,df$Patient_id)})
position = lapply(position,function(x){c(x,rep(NA,1-length(x)))})
position = unlist(position)
os_data$Position_matched = position

```

Now, we are creating a matrix `m` that we will fill with matching values

```
m = data.frame(matrix(ncol = length(df),nrow = nrow(os_data)))
colnames(m) = colnames(df)

for(i in 1:length(position))
{
  if(is.na(position[i])){}
  else{
    m[i,] = df[position[i],]
  }
}

```

And now, we can bind both `os_data` and `m`:

```
New_df = cbind(os_data,m)

```

And here the final output:

```
> New_df
   Patient_ids OS_time Position_matched ETV7 PLAUR SLC11A1                   Patient_id
1 TCGA-IG-A4QS    3.93               NA   NA    NA      NA                         <NA>
2 TCGA-JY-A6FB   61.23                1 1.79  6.24    3.70 TCGA-JY-A6FB-01A-11R-A336-31
3 TCGA-JY-A938   35.33                2 5.38  7.05    3.12 TCGA-JY-A938-01A-11R-A37I-31
4 TCGA-JY-A939   22.00               NA   NA    NA      NA                         <NA>
5 TCGA-JY-A93C   23.50               NA   NA    NA      NA                         <NA>
6 TCGA-JY-A93E   25.56               NA   NA    NA      NA                         <NA>
7 TCGA-L5-A4OJ   21.30               NA   NA    NA      NA                         <NA>
8 TCGA-L5-A4OS   59.40                3 5.00  5.00    5.00 TCGA-L5-A4OS-01A-11R-A212-31

```

**Warning:** This script should work for some simple case like yours but take extra cautious to the `grep` output as the short `patient_id` term of the `os_data` could eventually generate several hits instead of one.
> 2 votes

---
Tags: r
---


# Question
Title: How to extract hypervariable region of mitochondrial dna from its fasta file?
Recently I've done some analysis on human mitochondrial DNA and now I want to run the same analysis on only hypervariable regions of mitochondrial DNA to see whether the changes that I have observed are linked to only hypervariable regions. To do that, I'll need to make new fasta files from the mitochondrial DNA entries that I have and I want to know how to extract those. Any leads will be appreciated.

Link to one of the reference fasta files that I have

# Answer
> 2 votes
You can use Mitomap to get the coordinates of the hypervariable regions in the rCRS mitochondrial sequence (the standard MT sequence in the field). I checked the Genome Loci section of Mitomap and saw that it lists three:

```
Map Locus   Starting    Ending  Shorthand   Description
MT-HV1      16024       16383   CR:HVS1/HV1 Hypervariable segment 1
                                            [classic:16024-16365]
MT-HV2      57          372     CR:HVS2     Hypervariable segment 2
                                            [classic:73-340 ]
MT-HV3      438         574     CR:HVS3     Hypervariable segment 3

```

That gives you the coordinates on the rCRS sequence. What you can do now, is align all of your sequences and the rCRS together and see if you can infer where the hypervariable regions fall in your target sequences. Since, presumably, the rest of the sequence will be more conserved since it isn't hypervariable, you should be able to identify the syntenic blocks and locate the hypervariable regions in your sequences. I would expect to see a very high level of sequence identity in the non-hypervariable regions across all of your sequences (since they're all human mitochondria), so you should be able to find the hypervariable ones as the regions of low sequence similarity which also coincide with the known coordinates of the rCRS.

---
Tags: fasta, phylogenetics, database, phylogeny, genomics
---


# Question
Title: How to Installing & Using RDP Classifier through bioconda?
I am interested in using RDP classifier and RDPTools at the command line. In following the docs from RDP, I ran into this issue during the build, after running `make`:

```
BUILD FAILED
/home/denir/bin/RDPTools/Clustering/nbproject/build-impl.xml:896: The following error occurred while executing this line:
/home/denir/bin/RDPTools/Clustering/nbproject/build-impl.xml:1470: The following error occurred while executing this line:
/home/denir/bin/RDPTools/AlignmentTools/nbproject/build-impl.xml:896: The following error occurred while executing this line:
/home/denir/bin/RDPTools/AlignmentTools/nbproject/build-impl.xml:1414: The following error occurred while executing this line:
/home/denir/bin/RDPTools/ReadSeq/nbproject/build-impl.xml:920: The following error occurred while executing this line:
/home/denir/bin/RDPTools/ReadSeq/nbproject/build-impl.xml:300: Compile failed; see the compiler error output for details.

Total time: 1 second
Makefile:15: recipe for target 'Clustering/dist/Clustering.jar' failed
make: *** [Clustering/dist/Clustering.jar] Error 1

```

As one user suggested, I have tried downgrading from ant 1.10 to ant 1.9, this didn't work for me.

I would rather not downgrade to JDK 7 as another user suggested if possible so I am looking for alternatives. One user was able to install rdp classifier with bioconda. I was able to do that OK. However, I'm not sure where to go from here? The tutorial provided with commands like:

```
java -Xmx2g -jar ~/RDPTools/Clustering.jar derep -m '#=GC_RF' -o derep.fa all_seqs.ids all_seqs.samples alignment/*.fasta

```

Depend on the build created using the install method recommended by RDP.

Does anyone else use rdp classifier installed via conda? Is there a quickstart guide out there that I could refer to? Thanks so much for your help!

# Answer
> 1 votes
The benefit of using the conda package is that you no longer need to know or care exactly where the jar files are (that many java-based packages require this is TERRIBLE design). The example `java -jar ...` command in your post becomes:

```
Clustering derep -m '#=GC_RF' -o derep.fa all_seqs.ids all_seqs.samples alignment/*.fasta

```

in the bioconda version. In bioconda packages, there are wrapper (usually python) scripts meant to provide somewhat friendlier interfaces to java programs that otherwise lack executables. In this case, RDP Classfier has the following executables in bioconda:

* AbundanceStats
* AlignmentTools
* Clustering
* FrameBot
* KmerFilter
* ProbeMatch
* ReadSeq
* SeqFilters
* SequenceMatch
* hmmgs
* classifier

If you really want to know, the .jar file is under `$PREFIX/share/rdptools-$PKG_VERSION-$PKG_BUILDNUM/`, where `$PREFIX` is where ever the environment is and the other two variables indicate the package version and build number (at the moment, `rdptools-2.0.2-1` would be the likely folder name).

---
Tags: bioconda
---


# Question
Title: Can foldx position-scan produce different result with different lists of mutations?
I am running foldx position scan for GFPmut2 (PDB: 6GO9) and I use a config file with the following template:

```
command=PositionScan
pdb=6GO9.pdb
pdb-dir= <pdb location>
output-dir= <output directory>
output-file= <file name>
positions= <list of positions, comma separated: QA204R,KA107R...>

```

I had created three different config files with different position list.

Some of the positions were present in more than one files. Foldx produced different estimates of ΔΔG for these positions from these different config files. Example:

```
File1-LYSA41K   0
File1:LYSA41R   -0.353758
--
File2-LYSA41K   0
File2:LYSA41R   -0.0356187
--
File3-LYSA41K   0
File3:LYSA41R   -0.0031362

```

Another example

```
File1-LYSA52K   0
File1:LYSA52R   0.0057333
--
File2-LYSA52K   0
File2:LYSA52R   -0.00809613
--
File3-LYSA52K   0
File3:LYSA52R   0.0077189

```

I can certainly run foldx with a union of the three list of positions but I don't understand how this difference can even occur. Shouldn't each position scan be independent?

# Answer
**Minute values.** So just a reminder that the ∆∆G units is a kcal/mol. Water bounces around and collides with different energies following a Maxwell-Boltzmann distribution, where the average collision at 25°C is 0.59 kcal/mol (thermal energy constant). A hydrogen bond reduced the ∆G by 1–3 kcal/mol. Other interactions such as π-π, π-S are on the 1–2 mark. So the difference you are seeing is totally negligible.

**Heuristic.** The reason why you see the differences is that the force field calculations are heuristic and are not testing all perturbations —it is not only sampling a limited amount of conformers of side chains (kind of like PyMOL would be if it were FF based), but neighbour perturbations. So sometimes you hit a wee energy minima. Hence why in some protocols replicates are required —not mutation scanning mind you.

**Bonus: Protein pI.** Lastly, a caveat about pI of protein, you probably already know this, but best be safe. Your scores are all low, so I am guessing these surface mutations? Improving thermal stability is great for increasing the T\_M measured by DSC or catalysis if you are scanning a transition state as done for theozyme design, but it does not mean you get better solubility. In this case, you are introducing positive charges, which will alter the pI of the protein upwards and most protein have a pI less than 7 and the closer it is to pH 7 the more they crash out. If you are designing a high pI GFP, `Rosetta supercharge` does FF calculations to alter the pI of protein, which might be worth a gander.
> 1 votes

---
Tags: protein-structure, pdb
---


# Question
Title: Seurat VlnPlot presenting expression of multiple genes in a single cluster
Seurat VlnPlots are most commonly used to visualize differences in any given gene expression across multiple clusters or cell types. For example:

```
VlnPlot(object = MouseCellAtlas, idents = c("T cell","Neutrophil","Erythroblast","Monocyte","Macrophage"), features = 'Fzd9')

```

But suppose I want to view the expression of several genes in just `idents = "T cell"`. I can do this to generate three separate plots:

```
VlnPlot(object = mca, idents = 'T cell', features = c('Fzd9','Ctnnb1','Apc'), combine = TRUE, ncol = 3)

```

But ideally I'd like this second plot to look just like the first plot, except with only a single tissue on the legend ("T cell") and multiple genes in the various colored violins. Unfortunately, I am not aware of any simple way to combine these violins into a single plot area.

This doesn't do what I need:

```
VlnPlot(object = mca, idents = 'T cell', features = c('Fzd9','Ctnnb1','Apc'), combine = TRUE, ncol = 3)

```

Is there built-in functionality in Seurat for generating a VlnPlot for multiple genes in a single tissue in a single plot? Or can anyone suggest a workaround?

# Answer
To do so one workaround it to have your data in "long format" and then use the column that holds the "gene names" as the `x` variable while plotting.

You can use `FetchData()` to extract data from a `Seurat` object. `VlnPlot`'s default is the `data` slot (of the active assay if using Seurat v3 I suppose). And you can specify which cells and genes to retrieve.

```
selected_cells <- names(panc8$celltype[panc8$celltype == "gamma"])
data <- FetchData(panc8,
                    vars = c('FZD9','CTNNB1','APC'),
                    cells = selected_cells ,
                    slot = "data")
> head(data)
           FZD9   CTNNB1      APC
D101_5        0 0.000000 0.000000
D101_43       0 2.007853 1.001958
D101_93       0 0.000000 0.000000

```

`Melt()` would transform your data into "long format". By not specifying any arguments, we push for all of the info in the three variables to be gathered in two columns:

```
long_data <- melt(data)
No id variables; using all as measure variables
> head(long_data)
  variable value
1     FZD9     0
2     FZD9     0
3     FZD9     0

> tail(long_data)
     variable value
1870      APC     0
1871      APC     0
1872      APC     0

```

`ggplot2` is used to add "violin" and "jitter" layers. You can customize the output to look (exactly) like the `VlnPlot()` output.

```
ggplot(long_data,
       aes(x = variable, y = value)) +
  geom_violin() +
  geom_jitter(size = 0.1)

```

And here is the same graph generated by `VlnPlot()`:

There are no "violins" as the counts are almost entirely zeros. And see how different the ranges are in the y axes of the `VlnPlot`. So the "tweak" I have presented here would only work for genes that are expressed at similar levels / similar ranges.
> 2 votes

---
Tags: r, scrnaseq, seurat, ggplot2
---


# Question
Title: How to convert mzXML to mzML format (mass spectrometry files)?
How to convert mzXML to mzML format (mass spectrometry files) preferably under Linux and/or Python?

# Answer
`ProteoWizard` offers this functionality:

```
msconvert your_mzXML_file --mzML

```

`FileConverter` of `OpenMS` (I believe using the integrated ProteoWizard) also does this.
> 4 votes

---
Tags: python, linux
---


# Question
Title: Uploading files from local computer or web in rstudio cloud
I am using RStudio Cloud for running velocyto R package. So I have tried to import velocyto file from http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom

I tried to upload that like so

But saying

```
> SCG71 <- read.table("/cloud/project/SCG71.loom", header=TRUE, quote="\"")
Error in make.names(col.names, unique = TRUE) :
  invalid multibyte string at '<89>HDF'
In addition: There were 19 warnings (use warnings() to see them)

```

Can you help me in very initial of velocyto R?

# Answer
This seems to be an encoding problem (certain characters not being recognized), see this, this and also this. Probably the file `"/cloud/project/SCG71.loom"` is encoded with a format that is not included in your "locale".

Check the output of `Sys.getlocale()` and see if your input file's encoding is listed. If not, you will need to try something like `Sys.setlocale("LC_CTYPE", "en_US.UTF-8")` or `Sys.setlocale("LC_ALL", "English")` (see the third link for details).
> 1 votes

# Answer
**Velocyto.R**

I'm not sure that your issues is from RStudio cloud but more about how you are trying to open the file. The file that you are trying to import is not a `velocyto` file but a `.loom` file that can be open using the `velocyto.R` package.

I'm not very familiar with this kind of of `.loom` file but on the `velocyto.R` tutorial, it is specified to do the following:

```
library(velocyto.R)
ldat <- read.loom.matrices("SCG71.loom")

```

And you have the following file structures:

```
> str(ldat)
List of 3
 $ spliced  :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:2822201] 74 130 178 192 202 223 225 258 353 418 ...
  .. ..@ p       : int [1:6668] 0 464 786 1469 2502 3290 3760 4974 5388 5659 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:2822201] 1 1 1 1 1 1 1 1 1 4 ...
  .. ..@ factors : list()
 $ unspliced:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:528048] 61 100 339 352 377 445 465 528 617 758 ...
  .. ..@ p       : int [1:6668] 0 115 188 283 431 580 696 994 1127 1232 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:528048] 1 3 1 1 1 1 1 1 1 1 ...
  .. ..@ factors : list()
 $ ambiguous:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:47846] 11915 13649 13830 15638 18068 19585 23490 24345 8142 8767 ...
  .. ..@ p       : int [1:6668] 0 8 13 25 42 56 64 87 97 100 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:47846] 2 1 1 1 1 2 1 2 1 1 ...
  .. ..@ factors : list()

```

**Alternative 1 - LoomExperiment**

On bioconductor, you can also find `LoomExperiment`. On my computer, once I get it install, I was able to load your `SCG71.loom` file using:

```
library(LoomExperiment)
scle <- import("SCG71.loom",type = "LoomExperiment", rownames_attr = NULL, colnames_attr = NULL)

```

Then, you can find the output here:

```
> head(scle)
class: LoomExperiment
dim: 6 6667
metadata(4): chunks creation_date velocyto.__version__ velocyto.logic
assays(4): matrix ambiguous spliced unspliced
rownames: NULL
rowData names(6): Accession Chromosome ... Start Strand
colnames: NULL
colData names(1): CellID
rowGraphs(0): NULL
colGraphs(0): NULL

```

You can then load some data

```
> scle$CellID[c(1:5)]
[1] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG"
[4] "SCG71:AAACAAACAGCAGAAC" "SCG71:AAACAAACTTTGGGAG"

```

After that, I don't know what to do with that but I think your issue of importing a datafile is because you don't have the right package to open it.

**Alternative - using LoomR**

Alternatively, you can also install the package `loomR` and then enter:

```
library(loomR)
lfile = connect(filename = "SCG71.loom", mode = "r+")

```

And you get the output:

```
> lfile
Class: loom
Filename: /***my_directory***/SCG71.loom
Access type: H5F_ACC_RDWR
Attributes: creation_date, chunks, velocyto.__version__, velocyto.logic, version
Listing:
       name    obj_type dataset.dims dataset.type_class
  col_attrs   H5I_GROUP         <NA>               <NA>
 col_graphs   H5I_GROUP         <NA>               <NA>
     layers   H5I_GROUP         <NA>               <NA>
     matrix H5I_DATASET 6667 x 24421          H5T_FLOAT
  row_attrs   H5I_GROUP         <NA>               <NA>
 row_graphs   H5I_GROUP         <NA>               <NA>

```

Then, for manipulating these files, you should refer to their tutorial (it's beyond my skills and the scope of your question).

Hope it helps you to figure the solution of your question
> 1 votes

---
Tags: r, single-cell, rstudio, cloud
---


# Question
Title: Why are my genes filtered for Gene Ontology term enrichment?
When I enter a list of gene names into common gene ontology (GO) term enricher (e.g. MsigDB or ToppGene) only a small fraction of genes is actually used for the enrichment:

**I am struggling to understand why precisely only such a small number of genes is mapped here.**

My dataset is mapped with the GENCODE annotation. What annotation is the base of gene ontology? (The picture mentions entrez genes)

# Answer
What is "click here for details" saying? That link should contain an explanation of what happened.

GSEA usually requires `entrez_id` to run. If your list contains `gene_symbols`, then before enrichment, you need to convert the genes names. This can and is (always) tricky as genes are known with different names. For example, by googling I couldn't find any `AASS3` gene, but I found `AASS`, which is also named `LKRSDH`,`LORSDH`, and `LKR/SDH`. Depending on the software you are using for conversion, and the genome version used some genes will be missed. I usually run (in R) a combination of `geneSynonym` and `biomaRt` to fill the gaps in the lists.

Even better, would be to use `ensembl_id` or even better `entrez_id` directly
> 1 votes

---
Tags: gene-ontology, functional-annotation, go-enrichment, sequence-annotation, gencode
---


# Question
Title: How to convert a Pileup file to VCF format with Hg19 alignment
I would like to know how to convert a pileup file for an adna y chromosome to vcf with Hg19 alignment.

# Answer
> 2 votes
## Short answer

You'll have better luck searching for the term ***variant calling*** than with the term ***converting***. A Google search for "variant calling" gives several beginner-level tutorials on the topic. Commonly used tools for variant calling include samtools, GATK, and FreeBayes.

## Long answer

Converting data files from one format to another is a common task in bioinformatics. But that is not the task you want to perform here—the term ***convert*** is misleading.

A pileup file contains reads aligned to a reference sequence. A VCF file contains a description of putative or validated genetic variant. The process of creating a VCF from alignments is commonly called *variant discovery*, *variant prediction*, or *variant calling*. This is a broad and active area of research interest in bioinformatics, with many available algorithms and tools, each with a different set of parameters and assumptions.

This forum is not well designed to give you a good introduction to this topic. We are much better equipped to answer specific questions you have about how a particular tool or algorithm works, or why you are observing a particular unexpected result, and so on...assuming you have put in some time to research the question yourself.

---
Tags: samtools
---


# Question
Title: Proof of Breakpoint Reversal Sorting Approximation Algorithm
I'm kind of new to bioinformatics and trying to self-study. I'm reading a bioinformatics book: *An Introduction to Bioinformatics Algorithms* and ran into some problems about understanding the proof of the Breakpoint Reversal Sort algorithm in chapter 5.4.

Specifically, the book states a theorem:

> If a permutation $\pi$ contains a decreasing strip, then there is a reversal $\rho$ that decreases the number of breakpoints in $\pi$, that is, breakpoint ($\rho(\pi)$) \< breakpoint($\pi$).

This theorem guarantees that the algorithm will terminate and lead to no breakpoint. The book later explains the proof for the theorem:

> Let $k$ be the smallest element of decreasing strip, so element $k-1$ would be the ending of an increasing strip. Therefore element $k$ and $k-1$ corresponds to a breakpoint. Then reversing the segment between $k$ and $k-1$ will decrease the number of breakpoints.

I'm quite confused about this proof and don't really understand what the authors mean by reversing the segment between $k$ and $k-1$. Can someone please explain the proof for me?

Here's the lecture slide about reversal sorting that based on the book http://csbio.unc.edu/mcmillan/Media/Lecture07Spring2015.pdf.

# Answer
> 1 votes
Denote the strips in your permutation as $\langle S\_1, S\_2, \cdots, S\_n \rangle$. Suppose $k-1$ is in segment $S\_i$, and $k$ is in segment $S\_j$, for some $i, j \in \[n\], i \neq j$.

Recall by definition of $k$ that $k-1$ must be the maximum and therefore rightmost element in an increasing strip, and $k$ the minimum and therefore also rightmost in a decreasing strip. Consider the two cases:

* $i \< j$. The increasing strip comes before the decreasing. $k-1$ is on the 'right' end of $S\_i$, but $k$ is on the 'right' end of $S\_j$ also. In order for $k-1$ and $k$ to be adjacent, reverse the segments $S\_{i+1}, \cdots, S\_j$.
* $j \< i$. The decreasing strip comes before the increasing. $k$ is on the 'right' end of $S\_j$, but $k-1$ is on the 'right' end of $S\_i$ also. Again, we reverse $S\_{i+1}, \cdots, S\_j$.

So, "reversing the segment between $k$ and $k-1$" means to reverse all the strips in-between the two strips containing $k$ and $k-1$, excluding the 'leftmost' strip, but including the 'rightmost' strip.

---
Tags: sequence-alignment, genomics
---


# Question
Title: Trim Seurat object to contain expression info only for selected genes
I'm working with some large Seurat objects (MOCA, MCA, Tabula Muris) studying gene coexpression, and I'm running into memory issues.

Is it possible to remove all genes in a Seurat object that are not in a specified gene array? Ideally I'd like to be working with a small R matrix with say, 250k rows (cells) and 10-50 columns (gene expression data) which will enable fast processing and comprehensive analysis using other packages.

# Answer
> 1 votes
Is your memory usage low enough that you can initially read everything in and then filter later? If so you could read the Seurat object in like normal and then use `subset` against a vector of your gene names.

```
subsample <- subset(seurat_obj, features = my_genes)

```

# Answer
> 2 votes
A different approach if you are using `Seurat3`, is `DietSeurat()`. It allows you to diet the object by removing the components that you don't need. For example you can keep the normalised/scaled matrix and remove the raw counts. This approach could reduce space and memory usage, while keeping all your genes in place.

---
Tags: r, scrnaseq, seurat
---


title: Loss and gain of region in ATAC seq

text:
# Question
Title: Loss and gain of region in ATAC seq
I would like to calculate what is the number of region which is lost or gained between two condition,such as Stem cell vs Leukemic stemcell, atac seq data.

How do I get the information after peak calling such as what tools should i be using ?

As shown in the paper figure 1c I would like to do the same. How do i calculate the gain and loss of peak relative to a condition paper

# Answer
What you are looking for is "Differential accessibility analysis", you are using ATAC-seq, hence "accessibility". You will need to "call" and then "normalize" "count"(\*) peaks, after these, you can use `DESeq2` or `edgeR` for the "differential" part.

Here is a nice Nextflow workflow that addresses all of the steps of a differential accessibility analysis of ATAC-Seq data.

* `featureCounts` can be used the identify peaks, see this and this.
> 3 votes

---
Tags: atac-seq
---


title: How I can run this code on my files?

text:
# Question
Title: How I can run this code on my files?
I am annotating some .txt files by Annovar software by this code

```
nnovar]$ module load annovar/2016Feb01

[cyan01 annovar]$ table_annovar.pl

```

But I really got confused to do that by `awk` for all of my files in the path; Can somebody help me?

Many thanks

> EDITED

I found this but not working

```
module load annovar/20

for file in *.vcf; do
    table_annovar.pl example/a.txt humandb/ -buildver hg19 -out file -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f "$f" > "${f%.*}.txt"
done;

```

# Answer
`$f` in your command is not pointing anywhere since your for loop was defined as `for file in *.vcf;`. You should use `$file` instead of `$f`.

EDIT AS THE OP HAS EDITED THE QUESTION:

First of all, please copy and paste the error messages you are getting when seeking answers from the community, most of the time the error messages are clear enough to point to the right direction.

In you case, the problem is that you are trying to use a variable that you did not declare in the first place. It is explained above in my original answer but here I go again:

With this bit `for file in *.vcf;`, you are creating a variable `file` that you will using throughout the iteration. But then within the iteration, the next line (`... -operation gx,r,f,f,f "$f" > ...`, you try to use another variable `file`, that has not been declared yet (as far as I can see from your code bit).
> 4 votes

# Answer
The basic outline for loops is

`for sample in *; do mycommand "$sample" ; done`

You might need some way of incorporating the name of `$sample` into the output file name, so you don't keep overwriting the same output file over and over again (or do `mkdir` with each sample name and write everything for each sample in its own folder).
> 1 votes

---
Tags: bash, linux, shell, looper
---


# Question
Title: Fastq-dump script download X spots or all
Im trying to write a script where an optional input of -X flag can be used, or if that info is not available download all reads.

my script as follow:

```
# usage : bash dump_qc.sh <run id> <number of reads>

#run id variable.
SRR=$1

#maximum reads variable
X=$2

fastq-dump --split-3 -X ${X} ${SRR}

```

So, if I run the script

```
bash dump_qc.sh SRR1234 100

```

Will download 100 reads from the specified SRR file.

What I would like to do is to make the second argument optional, if not specified download the whole data set!

# Answer
> 1 votes
You don't need to copy your arguments into variables, you can just use `$1` directly. Anyway, this should be enough:

```
#!/bin/bash
# usage : bash dump_qc.sh <run id> <number of reads>

if [ -z $1 ]; then
  echo "ERROR: at least one argument need to be given!" >&2
  exit 1
fi

if [ ! -z $2 ]; then
    x="-X '$2'"
fi

fastq-dump --split-3 "$x" "$1"

```

---
Tags: fastq, bash, shell, sratoolkit
---


# Question
Title: Is there a Python package to convert InChi to molecular structures?
I am looking for a python package that can convert InChi keys (e.g.SGNXVBOIDPPRJJ-UHFFFAOYSA-N) to molecular structure or SMILES strings. Can I do that with RDKit?

The key comes from: https://www.metabolomicsworkbench.org/databases/refmet/refmet\_latest.txt

# Answer
> 3 votes
Yes, RDKit can be used, however, if you installed it with `conda` it will not work out of the box for inChi key fetching. You can either spend some time installing the missing bit or use something else. It is not like Mol2 support in RDKit, which shouldn't be used, it is just not worth the bother.

I would suggest PubChem's RESTFul API.

Say

```
 import requests
 ethanol = 'LFQSCWFLJHTTHZ-UHFFFAOYSA-N'
 r = requests.get(f'https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/inchikey/{ethanol}/property/CanonicalSMILES/JSON').json()
 smiles = r['PropertyTable']['Properties'][0]['CanonicalSMILES']

```

# Answer
> 2 votes
Slightly modified Matteos answer:

```
import requests
def get_smiles_from_inchikey(inchikey):
    r = requests.get(f'https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/inchikey/{inchikey}/property/CanonicalSMILES/JSON').json()
    return r['PropertyTable']['Properties'][0]['CanonicalSMILES']

inchikey = 'SGNXVBOIDPPRJJ-UHFFFAOYSA-N'
get_smiles_from_inchikey(inchikey)

```

---
Tags: structural-biology
---


# Question
Title: Putting cells in Seurat dimplot selectively
I am going to adjust Seurat dimplot in a way avoiding some cells so both my dimplot and heatmap look nice

If you please consider this picture, you would see some cells are far from the clusters so I want to avoid them in dimplot and of course for heatmap (coming from finding markers).

\[!\[enter image description here\]\[1\]\]\[1\]

Can you suggest me a piece of code doing that? Thanks

# Answer
> 1 votes
If you're using a GUI you could select the cells interactively:

```
plot <- DimPlot(seurat_obj, reduction = "umap")

```

Then select the cells by clicking around them

```
select.cells <- CellSelector(plot = plot)
Idents(seurat_obj, cells = select.cells) <- "SubCells"

```

and subset based on these cells

```
sub_cells <- WhichCells(seurat_obj, idents = "SubCells")
sub_obj <- subset(seurat_obj, cells = sub_cells)

```

You could then use subobject in your downstream dimplot.

See here for some of the interactive options. I guess you could outright subset on UMAP coordinates, but that may not always be easy if your boundary isn't linear.

---
Tags: r, seurat
---