# Question
Title: Determining position of side chain hydrogen in glycine residues from coordinates of backbone atoms
I have a dataset of pdb-format protein structure decoys/models generated by tertiary structure prediction programs. Many of these are missing hydrogen atoms.

Specifically, I'd like to determine the rough location of side-chain hydrogens in glycine residues, given the CA, C, N, and O coordinates. Is this possible? I don't need too much precision - I just need the general direction relative to the rest of the residue. More specifically, I want a vector originating at the CA that points roughly in the direction of the side chain.

I've looked at software for generating hydrogen atoms from pdb structures, however, I need to know which ones are glycine side chains, and was wondering if there might be a simpler way to predict hydrogen locations given coordinates of backbone atoms.

# Answer
> 3 votes
Glycine has a single hydrogen atom as its side chain:

All the six bond angles with the CA atom in the middle are about 109° (C-CA-N, C-CA-HA3, C-CA-HA2, N-CA-HA3, N-CA-HA2 and HA3-CA-HA2 using the CCD naming convention). This defines the rough direction of the hydrogen atoms.

But which of the two hydrogens is the side chain?

If it is a left-handed protein, I'd assume that it's the hydrogen pointing toward the screen, similarly to the side chain in L-alanine:

The images are taken from WikiMedia Commons: glycine and L-alanine.

---
Tags: proteins, protein-structure, pdb, 3d-structure
---


title: Issues installing Dendropy

text:
# Question
Title: Issues installing Dendropy
I want to perform phylogenetic analysis using Dendropy. I have installed dendropy using command `sudo pip install -U dendropy` but when I import it in python3.5.1 it gives followinng error:

```
 > python3
   Python 3.5.1 |Anaconda 2.5.0 (64-bit)| (default, Dec  7 2015, 11:16:01)
   [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] on linux
   Type "help", "copyright", "credits" or "license" for more information.
 >>> import dendropy
 Traceback (most recent call last):
 File "<stdin>", line 1, in <module>
 ImportError: No module named 'dendropy'

```

I again installed dendropy using the above given command as follows:

```
sudo pip3 install -U dendropy
Requirement already up-to-date: dendropy in /usr/local/lib/python3.5/dist-packages (4.4.0)
Requirement not upgraded as not directly required: setuptools in /usr/local/lib/python3.5/dist-packages (from dendropy) (39.2.0)

```

After doing this again the error remains the same and I am still not able to import dendropy, any solution to this problem?

# Answer
> 3 votes
You should **never** install python packages as root or with `sudo`.

You're already using conda, so `conda install -c bioconda -c conda-forge dendropy`.

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


# Question
Title: Problems in creating desired phylogenetic tree with ggtree
I am working on haplotype data and want to make a tree out of haplogroups using **ggtree**. I have following data in **newick format**.

```
(B(B2b_M112),C(C_RPS4Y711,C1a_V20,C2_M217,C2b1a_F3985,C2C1_Z1338),E(E1b1a_V38,E1b1b_M215,E2_M75),G(G_M201,G1_M342,G2a1_Z6552,G2a2a_PF3147,G2a2b2b_PF3359,G2a2b1_M406,G2a2b2a1a_U1,G2a2b2a1b_L497,G2a2b2a1c_CTS342,G2a_L30,G2b_M3155,G2b_M283),H(H1_M69,H1a_M82),I(I_M170,I1_M253,I2a1_Isles,I2a1_S21825,I2a1a_M26,I2a1b3_L621,I2a2a_M223,I2a2b_L38,I2c2_Y16419),J(J_M304,J1a_PF7257,J1a2a1a2_P58,J1a2a2_PF7264,J1a3_Z1828,J2_M172,J2a_L25,J2a1_L26,J2a1_M319,J2a1_Z1846,J2a1_Z387,J2a1_Z467,J2a1_Z6046,J2a1_Z6063,J2a1_Z6065,J2a1_Z7671,J2a1_Z7700,J2a2_PF5008, J2b_M241,J2b1_M205,J2b2_M241,J2b2_Z2456),K(K_M9),L(L_M20,L1_M22,L1a_M27,L1b_M317,L1c_M357),N(N3,N1a1_M46,N1a2_L666),O(O1a1_B384,O1b1_F2320,O1b2_P49,O2 IMS_JST0213554,O2a1_F51,O_O2a2_F525,O3_M122),Q(Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),R(R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),T(T_M184,T_L131,T_PF5633));

```

I am using following R code to plot a tree from this data,

```
> library("ggtree")
> tree <- read.tree(text = "(B(B2b_M112),C(C_RPS4Y711,C1a_V20,C2_M217,C2b1a_F3985,C2C1_Z1338),E(E1b1a_V38,E1b1b_M215,E2_M75),G(G_M201,G1_M342,G2a1_Z6552,G2a2a_PF3147,G2a2b2b_PF3359,G2a2b1_M406,G2a2b2a1a_U1,G2a2b2a1b_L497,G2a2b2a1c_CTS342,G2a_L30,G2b_M3155,G2b_M283),H(H1_M69,H1a_M82),I(I_M170,I1_M253,I2a1_Isles,I2a1_S21825,I2a1a_M26,I2a1b3_L621,I2a2a_M223,I2a2b_L38,I2c2_Y16419),J(J_M304,J1a_PF7257,J1a2a1a2_P58,J1a2a2_PF7264,J1a3_Z1828,J2_M172,J2a_L25,J2a1_L26,J2a1_M319,J2a1_Z1846,J2a1_Z387,J2a1_Z467,J2a1_Z6046,J2a1_Z6063,J2a1_Z6065,J2a1_Z7671,J2a1_Z7700,J2a2_PF5008, J2b_M241,J2b1_M205,J2b2_M241,J2b2_Z2456),K(K_M9),L(L_M20,L1_M22,L1a_M27,L1b_M317,L1c_M357),N(N3,N1a1_M46,N1a2_L666),O(O1a1_B384,O1b1_F2320,O1b2_P49,O2 IMS_JST0213554,O2a1_F51,O_O2a2_F525,O3_M122),Q(Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),R(R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),T(T_M184,T_L131,T_PF5633));")
> p<-ggtree(tree)
> p+geom_tiplab()

```

And the displayed tree looks like this

I have following problems in this tree,

1. It gives a low quality image
2. It does not show the haplotype taxa names clearly
3. I want to arrange haplogroups i.e nodes in the desired according to the tree given below,

i.e B haplogroups with C haplotypes and with C, E named haplotypes etc.

I have tried to solve these problems but failed to find any solution. Is there a solution?

# Answer
> 1 votes
The above given problems were solved using manual tree construction through MEGA. MEGA allows manual construction/manipulation of trees. The newick format obtained after manually constructing the desired tree using MEGA is:

```
((((((((((((T_M184,T_L131,T_PF5633),T),((R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),R)),((Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),Q)),(((((((O1a1_B384,O1b1_F2320),O1b2_P49),O2IMS_JST0213554),O2a1_F51),O_O2a2_F525),O3_M122),O)),(((N3,N1a1_M46),N1a2_L666),N)),(((((L_M20,L1_M22),L1a_M27),L1b_M317),L1c_M357),L)),(K,K_M9)),((((((((((((((((((((((J_M304,J1a_PF7257),J1a2a1a2_P58),J1a2a2_PF7264),J1a3_Z1828),J2_M172),J2a_L25),J2a1_L26),J2a1_M319),J2a1_Z1846),J2a1_Z387),J2a1_Z467),J2a1_Z6046),J2a1_Z6063),J2a1_Z6065),J2a1_Z7671),J2a1_Z7700),J2a2_PF5008),J2b_M241),J2b1_M205),J2b2_M241),J2b2_Z2456),J)),(((((((((((I_M170,I1_M253),I2a1_Isles),I2a1_S21825),I2a1a_M26),I2a1b3_L621),I2a2a_M223),I2a2b_L38),I2c2_Y16419),I),((H1_M69,H1a_M82),H)),((((((((((((G_M201,G1_M342),G2a1_Z6552),G2a2a_PF3147),G2a2b2b_PF3359),G2a2b1_M406),G2a2b2a1a_U1),G2a2b2a1b_L497),G2a2b2a1c_CTS342),G2a_L30),G2b_M3155),G2b_M283),G))),((E1b1a_V38,E1b1b_M215),(E2_M75,E))),(((((((C_RPS4Y711,C1a_V20),C2_M217),C2b1a_F3985),C2C1_Z1338),C),B2b_M112),B));

```

---
Tags: r, bioconductor, phylogenetics, phylogeny
---


title: Macs2 peak calling?

text:
# Question
Title: Macs2 peak calling?
I have paired end ChIP-seq data with 101 bp and 2 biological replicates for each one. I have done peak calling with macs2 but I have some questions about it.

I also faced with an warning:

> WARNING @ Thu, 07 Jun 2018 17:06:05: #2 Since the d (197) calculated from paired-peaks are smaller than 2\*tag length, it may be influenced by unknown sequencing problem!
> WARNING @ Thu, 07 Jun 2018 17:06:05: #2 You may need to consider one of the other alternative d(s): 197
> WARNING @ Thu, 07 Jun 2018 17:06:05: #2 You can restart the process with --nomodel --extsize XXX with your choice or an arbitrary number. Nontheless, MACS will continute computing.

1. I have added `--nomodel --extsize 197 ; --nomodel --extsize 147` and `--nomodel --extsize 202` (separately to macs2 command) and got the results without any warning? Which one is more correct?
2. Are broad peaks extended of narrow peaks? If I apply intersect between them I should expect find 100% overlap between narrow peaks and broad ones?
3. Which kind of peak (narrow/broad) is proper for H3k27ac, H3k4me1, H3k4me3, H3k27me3 study?
4. If there is no control group for using as background, can I use default parameters?

# Answer
> 4 votes
1. Your original command without `--nomodel --extsize ...` is probably the most accurate. This warning stems from a time when reads were much much shorter and likely never made that much sense to begin with.
2. Broad peak calling in MACS2 basically works by finding a bunch of nearish narrow peaks and merging them. If you have really broad signals then use something like `histoneHMM` instead.
3. H3K27ac, H3K4me1 and H3K4me3 are narrow. H3K27me3 is broad.
4. Yeah, but you'll need to spend a bit more time with the data to ensure that the peaks are reasonable (you'll have to tweak the parameters if not).

---
Tags: chip-seq
---


title: Blueprint replicates

text:
# Question
Title: Blueprint replicates
On Encode, it's pretty easy to find out replicate information about the experiments, however I am having some difficulty finding this information on BLUEPRINT DDC, at http://dcc.blueprint-epigenome.eu

Where can I find information which describes technical replicates and biological (isogenic) replicates?

Thanks in advance.

# Answer
> 1 votes
Blueprint does not have technical replicates, you will find that very few sequencing based experiments ever have technical replicates. As biological replicates: each individual is regarded as a replicate.

This makes sense for the sort of experimental design which blueprint was designed to answer questions such as "Do individuals from group A have higher or lower H3K4me3 than individuals from groups B" or "Does cell type A have higher accessibility than cell type B". In each of these cases the within group variability to be considered is encapsulated in the difference between individuals.

Overall technical variation is assumed to be similar in all samples or at least not correlated with the contrast of interest (this may well not be true in the second example design I gave) it is safe to directly measure it, as long we know its distribution.

As an aside, one of the reasons that ENCODE got away with only two reps is that their reps are effectively technical (repeated preparations from the same cell lines) and thus there is next to no biological variation (which is what our algorithms are set up to find).

# Answer
> 0 votes
Unfortunately this information isn't as clearly laid out in Blueprint as it is in ENCODE. You can view metadata for experiments from one of the many data viewers (IHEC Data Portal, the Blueprint DCC you mentioned, DeepBlue, or the EGA dataset itself), but these tend to not contain technical/biological replicate information like ENCODE does.

You often have to dive into the data from the original publication to get that information

---
Tags: database, data-retrieval
---


# Question
Title: Using R to import specific records from a PDB file
I'm using R and I'm trying to import specific records from a PDB file (e.g. `HETNAM`, `HETSYN`). Unfortunately, it seems the `read.pdb` function of either `biomaRt` or `Rpdb` aren't capable of it.

Any ideas?

Thanks

# Answer
> 4 votes
If your file is in PDB format as the one shown here in .pdb format then you can simply extract your desired fields using sed which is a programming language used to parse and transform text. As for thisstructurei.e 6C4A in .pdb format HETNAM and HETSYN can be extracted using sed as follows:

```
sed  -n -e '/HETNAM/ p' -e '/HETSYN/ p'  6C4A.pdb

HETNAM     EJA S-[(1Z)-2-CARBOXY-N-HYDROXYETHANIMIDOYL]-L-CYSTEINE
HETNAM     PYR PYRUVIC ACID
HETNAM      MG MAGNESIUM ION
HETNAM     3NP 3-NITROPROPANOIC ACID
HETNAM     PEG DI(HYDROXYETHYL)ETHER
HETNAM     EPE 4-(2-HYDROXYETHYL)-1-PIPERAZINE ETHANESULFONIC ACID
HETSYN     EPE HEPES

```

---
Tags: r, pdb
---


# Question
Title: Use of heterozygous SNPs in cancer research: why?
When reading about allelic fraction (AF) and SNPs in cancer research, they always mention the fact that they're using heterozygous SNPs (*informative* SNPs). Why is this? Why can't we use homozygous SNPs?

In this case, what is the reference base (from the reference genome)?

# Answer
I’m no longer working in tumour sequencing so I’m by no means an expert. But in a nutshell, the reason is that, as indicated, homozygous SNPs *aren’t informative*: if your allelic fraction is 100%, meaning you only observe a single nucleotide at a given position, we don’t know whether we’re dealing with a tumour related mutation. In fact, we probably don’t because the sample likely contains healthy as well as tumour cells.

By contrast, tumour-specific variants always have an AF \< 100%, in other words they exhibit some degree of heterozygosity.

In fact, in healthy tissues you’d expect to see one of three cases: your variant allele (= sequence base ≠ reference) has a frequency of either 0% or 100% (homozygous), or 50% (heterozygous), ignoring rare somatic mutations. With tumour samples, by contrast, the fraction will vary depending on the penetrance of the variant in the tumour, and on the fraction of healthy to tumour cells in the sample.
> 5 votes

# Answer
I'm not sure what I'm talking about, but I guess one could re-phrase Konrad's answer in terms of information theory.

The information content of a message (a packet of bits) is inversely proportional to the frequency of those bits. For example, if you randomly extract a word (the message) from an English book and the word is "the", you don't get much information about the book because "the" is very frequent and you see it everywhere. If, instead, the word happens to be the much rarer "bioinformatics" you start having a pretty good guess what the book is about.

About SNPs, homozygote reference genotypes have very little information because most of the genomic positions in most of the people are homozygote reference. Heterozygote and homozygote non-reference are rarer and thus more informative.

So, I don't think there is anything wrong with using homozygote SNPs. They are ignored because they are so common that their information content is very small.
> 2 votes

---
Tags: ngs, snp, cancer
---


# Question
Title: Disk space error while aligning reads using STAR
Hi am trying to align RNA-seq reads using STAR through following commands,

```
STAR --runMode alignReads --runThreadN 26  --outSAMtype BAM  SortedByCoordinate Unsorted --sjdbOverhang 100  --genomeDir wheat --readFilesIn G3_cleaned_R1.fastq G3_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript wheat.gff3
Jun 05 13:53:59 ..... started STAR run
Jun 05 13:53:59 ..... loading genome
Jun 05 14:07:12 ..... started mapping
Jun 05 18:42:47 ..... started sorting BAM

EXITING because of FATAL ERROR: number of bytes expected from the BAM bin   does not agree with the actual size on disk: 509186583   0   45

EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 648593301   0   47

Jun 05 18:42:47 ...... FATAL ERROR, exiting
Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 660412146   0   46
Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 470931698   0   44

Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 653155717   0   43

Jun 05 18:42:47 ...... FATAL ERROR, exiting
Segmentation fault (core dumped)

```

Although there is much free space on my system as shown below:

```
> df -Th
  Filesystem     Type      Size  Used Avail Use% Mounted on
  devtmpfs       devtmpfs   63G     0   63G   0% /dev
  tmpfs          tmpfs      63G  136K   63G   1% /dev/shm
  tmpfs          tmpfs      63G  1.9M   63G   1% /run
  tmpfs          tmpfs      63G     0   63G   0% /sys/fs/cgroup
  /dev/sda3      ext4      3.6T  1.4T  2.2T  39% /
  /dev/sda2      ext4      9.8G   57M  9.2G   1% /boot
  tmpfs          tmpfs      13G     0   13G   0% /run/user/480
  tmpfs          tmpfs      13G   16K   13G   1% /run/user/1000

```

How this problem can be solved?

# Answer
STAR versions 2.6.0b and 2.6.0a are unstable, as written Issues performing variant calling with GATK you are using version 2.6.0b. You should switch to version 2.6.0c which is stable. Also avoid to perform multiple alignments in the same directory.
> 4 votes

---
Tags: sam, rna-alignment, star
---


# Question
Title: What is "Possible polymerase run-on fragment"?
I used gffcompare recently to compare a newly assembled transcriptome and a published transcriptome. The class code "p" transcripts are "possible polymerase run on".

What does that mean? I tried very hard but still cannot find the answer.

# Answer
> 5 votes
The general idea behind a run-on fragment is that it's background noise. This derives from an open area of the genome that's next to an area that's actually being transcribed. Thus, all of the machinery for transcription is already present, so if the region adopts a bit more of an open conformation then that machinery can sporadically bind and start transcribing. You can get those either up- or down-stream of genes. A lot of things like this are just background transcriptional noise.

# Answer
> 6 votes
On most coding genes (with the exception of replication dependent histone genes), *Transcripts* (as opposed to transcrip**tion**) are terminated by cleavage polyadenylation sites, where the growing transcript is cleaved by a protein complex and the poly-A tail. However, PolII continues to transcribe until it is terminated by torpedo termination (see here). Normally this extra RNA after the polyA site is degraded by Xrn2, but if you capture any particular snapshot of the RNA in the cell at any one time there will be a small amount that corresponds to this run on transcription.

---
Tags: transcriptome
---


# Question
Title: What does an output of -I for all amino acids mean in a psi-blast pssm matrix file?
I have run psi-blast using the NR database, remotely, with one iteration, for several sequences to calculate an evolutionary profile (PSSM) for each of those sequences.

However, many of the PSSM files contain lines of `-I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I` rather than numerical values. What does it mean?

```
>     Last position-specific scoring matrix computed
>                 A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V
>         1 G  -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I
>         2 N     4  -1  -2  -2   0  -1  -1   0  -2  -1  -1  -1  -1  -2  -1   1   0  -3  -2   0
>         3 K    -2  -1   4   4  -3   0   1  -1   0  -3  -4   0  -3  -3  -2   0  -1  -4  -3  -3
>         4 E     0  -3  -3  -3   9  -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1
>         5 K    -2  -2   1   6  -3   0   2  -1  -1  -3  -4  -1  -3  -3  -1   0  -1  -4  -3  -3
>         6 A    -1   0   0   2  -4   2   5  -2   0  -3  -3   1  -2  -3  -1   0  -1  -3  -2  -2
>         7 D    -2  -3  -3  -3  -2  -3  -3  -3  -1   0   0  -3   0   6  -4  -2  -2   1   3  -1
>         8 R     0  -2   0  -1  -3  -2  -2   6  -2  -4  -4  -2  -3  -3  -2   0  -2  -2  -3  -3
>         9 Q    -2   0   1  -1  -3   0   0  -2   8  -3  -3  -1  -2  -1  -2  -1  -2  -2   2  -3
>        10 K    -1  -3  -3  -3  -1  -3  -3  -4  -3   4   2  -3   1   0  -3  -2  -1  -3  -1   3
>        11 V    -1   2   0  -1  -3   1   1  -2  -1  -3  -2   5  -1  -3  -1   0  -1  -3  -2  -2
>        12 V    -1  -2  -3  -4  -1  -2  -3  -4  -3   2   4  -2   2   0  -3  -2  -1  -2  -1   1
>        13 S    -1  -1  -2  -3  -1   0  -2  -3  -2   1   2  -1   5   0  -2  -1  -1  -1  -1   1
>        14 D    -2   0   6   1  -3   0   0   0   1  -3  -3   0  -2  -3  -2   1   0  -4  -2  -3
>        15 L    -1  -2  -2  -1  -3  -1  -1  -2  -2  -3  -3  -1  -2  -4   7  -1  -1  -4  -3  -2
>        16 V    -1   1   0   0  -3   5   2  -2   0  -3  -2   1   0  -3  -1   0  -1  -2  -1  -2
>        17 A    -1   5   0  -2  -3   1   0  -2   0  -3  -2   2  -1  -3  -2  -1  -1  -3  -2  -3
>        18 L     1  -1   1   0  -1   0   0   0  -1  -2  -2   0  -1  -2  -1   4   1  -3  -2  -2
>        19 E     0  -1   0  -1  -1  -1  -1  -2  -2  -1  -1  -1  -1  -2  -1   1   5  -2  -2   0
>        20 G     0  -3  -3  -3  -1  -2  -2  -3  -3   3   1  -2   1  -1  -2  -2   0  -3  -1   4
>        21 A    -3  -3  -4  -4  -2  -2  -3  -2  -2  -3  -2  -3  -1   1  -4  -3  -2  11   2  -3
>        22 L    -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
>        23 D    -2  -2  -2  -3  -2  -1  -2  -3   2  -1  -1  -2  -1   3  -3  -2  -2   2   7  -1
>        24 M    -1   0   0   1  -3   4   4  -2   0  -3  -3   1  -1  -3  -1   0  -1  -2  -2  -2
>        25 Y     0  -3  -3  -3   9  -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1
>        26 K    -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4
>        27 L    -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
>        28 D    -1  -2  -3  -3  -1  -2  -3  -4  -3   3   3  -3   2   0  -3  -2  -1  -2  -1   2
>        29 N  -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I
>        30 S    -2   1  -1  -1  -2  -1  -2  -2  -1  -1  -1  -1  -2   0  -1   0   0  -1  -1  -3
>        31 R     4   0  -2   1   4  -1  -3   0  -3   0  -3  -4   4  -1   0  -1  -3  -1   0  -4
>        32 Y    -3  -1  -3  -4  -3  -3  -2  -3  -1  -3  -1  -1  -3  -3  -3  -1  -1  -1  -3  -2

```

# Answer
> 1 votes
It means -inf. It has to do with how the PSSM score is calculated.

---
Tags: blast, sequence-homology, phylogenetics, pssm, psi-blast
---


# Question
Title: Test to determine if two genes/exons share the same evolutionary histories?
In classic phylogenetic inference one is usually given various orthologue sequences of a given gene across various species. Those sequences are then multiple aligned and used to construct a phylogenetic tree.

Say, I split the multiple alignment of this gene into two parts (part A and B). My question now is: What is the most appropriate test to determine if part A and B share the same phylogenetic history?

Having looked through literature, it appears this can be done by constructing phylogenetic trees separately for part A and B and test for incongruence using

* Kishino–Hasegawa (**KH test**);
* Shimodaira–Hasegawa (**SH-test**);
* **SOWH-test** etc.

as nicely summarized in Plantet et all, and Goldman et al. etc.

However the above test rely on constructed trees beforehand which might be incorrect due to inappropriate algorithms/parameters.

I was therefore wondering if there are more appropriate test maybe without having to construct trees at all?

# Answer
Since there hasn't been any other answers, I'd thought to share my approach in case someone is having the same issue.

Gauging from the comments of tendon and NatWH, it seems using the SH-test valid stategy to tackle this problem.

Having done some research, I could only find one functional software package, the R-package phangorn which implements both, the SH- as well as SOWH-test and thought to share my R-workflow:

**Reading in multiple-alignment (MLN). Here: amino-acids (AA)**

```
 library(phangorn)
 mln <- read.phyDat("my_mln.fasta", format = "fasta", type = "AA")

```

**Partitioning MLN into two parts (part A and B) given a predefined breakpoint (here 100)**

```
breakpoint <- 100
 # partA: 1:breakpoing
 ( mln_partA <- subset(mln, select = 1:breakpoint, site.pattern = FALSE) )
 # 58 sequences with 100 character and 45 different site patterns.
 # The states are a r n d c q e g h i l k m f p s t w y v

 ## partB: breakpoint:end
 ( mln_partB <- subset(mln, select = (breakpoint + 1):length(attr(mln, "index")), site.pattern = FALSE) )
 # 58 sequences with 1702 character and 720 different site patterns.

```

**Constructing tree(s) separately from part A and B using Neighbourhood Joining (NJ) method**

```
 ## Note, other tree building algorithms might be used
 ( treeNJ_partA <- NJ(dist.ml(mln_partA)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.
 ( treeNJ_partB <- NJ(dist.ml(mln_partB)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.

```

**Using topology of previously constructed trees and compute optimal likelihood based on partA data**

```
 ( treeNJ_partA_fit1 <- optim.pml(pml(treeNJ_partA, mln_partA)) )
 # optimize edge weights:  -1402.126 --> -1379.043
 #
 #  loglikelihood: -1379.041
 #
 # unconstrained loglikelihood: -302.6292
 ( treeNJ_partB_fit1 <- optim.pml(pml(treeNJ_partB, mln_partA)) )
 # optimize edge weights:  -1633.646 --> -1530.745
 #
 #  loglikelihood: -1530.682
 #
 # unconstrained loglikelihood: -302.6292

```

**Performing SH-test for for both trees**

```
 ( shPart1_t12 <- SH.test(treeNJ_partA_fit1, treeNJ_partB_fit1, data = mln_partA, B = 1e5) )
 #      Trees      ln L Diff ln L p-value
 # [1,]     1 -1379.041     0.000 0.48461
 # [2,]     2 -1530.682   151.641 0.00000

```

As can be seen in this example, the tree constructed for partB doesn't fit partA data as much as the tree constructed from partA.
> 2 votes

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


# Question
Title: How to interpret a PSSM matrix from PSI-Blast - what do the numbers represent?
PSI-Blast on a sequence returns a PSSM matrix file containing two matrices, one to the left, one the the right, for all amino acids:

```
A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V

```

What is the difference between these two matrices, how can I interpret them?

```
Last position-specific scoring matrix computed, weighted observed percentages rounded down, information per position, and relative weight of gapless real matches to pseudocounts
            A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V   A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V
    1 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    2 G     0  -3   0  -2  -3  -2  -2   7  -2  -4  -4  -2  -3  -4  -2   0  -2  -3  -4  -4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    3 G     0  -3   0  -2  -3  -2  -2   7  -2  -4  -4  -2  -3  -4  -2   0  -2  -3  -4  -4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    4 G     0  -3   0  -2  -3  -2  -2   7  -2  -4  -4  -2  -3  -4  -2   0  -2  -3  -4  -4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    5 L    -2  -3  -4  -4  -1  -2  -3  -4  -3   2   4  -3   2   0  -3  -3  -1  -2  -1   1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    6 C     0  -4  -3  -4  10  -3  -4  -3  -3  -1  -1  -4  -2  -3  -3  -1  -1  -3  -3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    7 Y    -2  -2  -2  -4  -3  -2  -2  -4   2  -2  -1  -2  -1   3  -3  -2  -2   3   8  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    8 C     0  -4  -3  -4  10  -3  -4  -3  -3  -1  -1  -4  -2  -3  -3  -1  -1  -3  -3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    9 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   10 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   11 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   12 F    -3  -3  -3  -4  -3  -4  -4  -4  -1   0   0  -4   0   7  -4  -3  -2   1   3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   13 C     0  -4  -3  -4  10  -3  -4  -3  -3  -1  -1  -4  -2  -3  -3  -1  -1  -3  -3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   14 V     0  -3  -3  -4  -1  -3  -3  -4  -4   3   1  -3   1  -1  -3  -2   0  -3  -1   4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   15 C     0  -4  -3  -4  10  -3  -4  -3  -3  -1  -1  -4  -2  -3  -3  -1  -1  -3  -3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   16 V     0  -3  -3  -4  -1  -3  -3  -4  -4   3   1  -3   1  -1  -3  -2   0  -3  -1   4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   17 G     0  -3   0  -2  -3  -2  -2   7  -2  -4  -4  -2  -3  -4  -2   0  -2  -3  -4  -4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   18 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00

```

Additionally, at the end of the file, there is additional information, what does it show exactly?

```
                      K         Lambda
Standard Ungapped    0.1628     0.3454
Standard Gapped      0.0410     0.2670
PSI Ungapped         0.1628     0.3454
PSI Gapped           0.0410     0.2670

```

# Answer
> 1 votes
If you read wiki carefully, you can see that PSSM is calculated in 3 steps. First frequency is calculated (how much times was the amino acid or nucleotide on that location in the motif), from that you can calculate probability (in the wiki example there were 10 sequences, so each freq. is divided by 10).

Then the log likelihood is calculated, which are the PSSM values. These are in your first matrix (rounded down). The second matrix shows how much the values are relative from your pseudocounts (I assume it was set to default = 0).

Lambda and kappa are estimated to calculate the normalized score (S') for HSP, if you never heard of this before I suggest you read the original psi-blast paper first which you can find here.

---
Tags: blast, phylogenetics, pssm, psi-blast
---


# Question
Title: Problem installing last build of Sniffles
I have followed the steps described here to install Sniffles on my Mac Book Pro running macOS High Sierra.

However, I am struggling with this step:

```
cmake -D CMAKE_C_COMPILER=/opt/local/bin/gcc-mp-4.7 -D CMAKE_CXX_COMPILER=/opt/local/bin/g++-mp-4.7 ..

```

I do not have an /opt folder, so I get the following output:

```
-- The C compiler identification is unknown
-- The CXX compiler identification is unknown
CMake Error at CMakeLists.txt:2 (project):
  The CMAKE_C_COMPILER:

    /opt/local/bin/gcc-mp-4.7

  is not a full path to an existing compiler tool.

  Tell CMake where to find the compiler by setting either the environment
  variable "CC" or the CMake cache entry CMAKE_C_COMPILER to the full path to
  the compiler, or to the compiler name if it is in the PATH.

CMake Error at CMakeLists.txt:2 (project):
  The CMAKE_CXX_COMPILER:

    /opt/local/bin/g++-mp-4.7

  is not a full path to an existing compiler tool.

  Tell CMake where to find the compiler by setting either the environment
  variable "CXX" or the CMake cache entry CMAKE_CXX_COMPILER to the full path
  to the compiler, or to the compiler name if it is in the PATH.

-- Configuring incomplete, errors occurred!
See also "/Users/xxx/Documents/xxx/xxx/Sniffles-master/build/CMakeFiles/CMakeOutput.log".
See also "/Users/xxx/Documents/xxx/xxx/Sniffles-master/build/CMakeFiles/CMakeError.log".

```

Unsurprisingly, the error file says that /opt/local/bin/gcc-mp-4.7 does not exist.

So I tried:

```
cmake -D CMAKE_C_COMPILER=/usr/bin/gcc -D CMAKE_CXX_COMPILER=/usr/bin/gcc ..

```

I get another long message but the most important bit is:

```
CMake Error at /usr/local/Cellar/cmake/3.8.1/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:137 (message):
  Could NOT find OpenMP (missing: OpenMP_C_FLAGS OpenMP_CXX_FLAGS)

```

I do not know what to do to solve this error.

I do have gcc here: /usr/bin/gcc and cmake version 3.8.1.

I have created an issue on the Sniffles GitHub repository.

### Update 13th June 2018

I installed the latest version of gcc:

```
brew install gcc
Warning: gcc 8.1.0 is already installed and up-to-date
To reinstall 8.1.0, run `brew reinstall gcc`

```

Then, I tried to make cmake use gcc and g++ in the /usr/local/bin folder.

```
cmake -D CMAKE_C_COMPILER=/usr/local/bin/gcc-8 -D CMAKE_CXX_COMPILER=/usr/local/bin/g++-8 ..

```

This recreated the same problem:

```
CMake Error at /usr/local/Cellar/cmake/3.8.1/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:137 (message):
  Could NOT find OpenMP (missing: OpenMP_C_FLAGS OpenMP_CXX_FLAGS)

```

Do I need to set the OpenMP\_C\*\_FLAGS to something?

# Answer
> 4 votes
You need to install GCC, or a newer version of Clang.

Your compiler — though the binary is called `gcc` — is actually Clang, and the Clang version that’s shipping with macOS/Xcode still does not support OpenMP.

The easiest course of action is probably to install GCC via Homebrew (`brew install gcc`) and then configure Sniffles to use Homebrew’s GCC (i.e. setting the `CMAKE_C*_COMPILER` paths to `/usr/local/bin/g{cc,++}-$version`).

---
Tags: software-installation
---


# Question
Title: Tumor sample and heterozygous SNVs: what is the genotype of the normal cells?
In the following picture, you can see that in normal cells (yellow) we have two copies of each gene, while in tumor cells (violet), we just have one copy due to a clonal mono-allelic deletion. I do not really care about the genotype of the two genes, I just need to know that they are heterozygous.

On the other hand, suppose we are talking about somatic heterozygous SNVs and therefore we do not consider deletions anymore. In this case I would like to know the genotype (AA, AB or BB) of tumor and normal cells. Based on my understanding, I expect the AF of the SNVs in the tumor cells to be roughly 0.5 (since I ask for heterozygous SNVs). But what about normal cells? What is their genotype? If I ask for heterozygous SNVs, does that imply that normal cells are homozygous for that alleles?

# Answer
Your initial information/image and actual question seem pretty unrelated to me, but I'll do what I can.

If you have a heterozygous SNV in a malignant population of cells, it really depends on when the mutation occurred for you to know the genotype of wildtype cells. It could have been there since birth (a **germline** mutation passed down from the parents) or have been acquired throughout the person's life and clonally expanded (**somatic** mutation).

Unfortunately, you have no way to differentiate between these two types of SNVs just by sequencing the tumor cells. If you want to know true **somatic** mutations that are found only in the tumor cells, you need to sequence a matched normal control from the same patient (preferably from the same cell type). This is why most large cancer sequencing studies include not only tumor samples, but normal cells from the same patient - see TCGA or ICGC for examples.
> 0 votes

---
Tags: ngs, cancer, snv
---


# Question
Title: Seurat CCA, AlignSubspace, Argument dist.method is only useful with multivariate timeseries
I am trying to combine two datasets to be able to compare the cells in them and using `Seurat` `CCA` for it:

https://satijalab.org/seurat/Seurat\_AlignmentTutorial.html

When I am running `AlignSubspace` function I am getting a warning for each single dimension that I selected it to be run on:

```
seurat_union <- AlignSubspace(object = seurat_union, reduction.type =
    "cca", grouping.var = "protocol", dims.align = 1:13)

```

So, I am getting 13 warnings:

```
> warnings()
Warning messages:
1: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
2: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
3: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
4: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
5: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
6: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
7: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
8: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
9: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
10: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
11: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
12: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
13: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries

```

What is this warning telling me? What should I pay attention to with my datasets?

# Answer
> 1 votes
I don't have enough reputation to comment so I will briefly answer here

This warning is fixed in Seurat 2.3.2. Without upgrading you can ignore this message. It doesn't mean anything for the actual analysis performed so your results won't be flawed.

---
Tags: seurat
---


# Question
Title: Running python script from terminal: umap issue
Some time ago I asked the question:

RunUMAP in Seurat not working: module 'umap' has no attribute 'UMAP'

It was solved by following the answer in it, however, it is still not working when I run a python script from terminal, giving the same error:

```
python python/umap.py -path
"/Users/nikitavlasenko/python_scripts/output/" -data
"/Users/nikitavlasenko/python_scripts/data_files/WT_SCID_CCA.csv"

```

> Traceback (most recent call last): File "python/umap.py", line 26, in result = umap.UMAP().fit\_transform(data.T) AttributeError: module 'umap' has no attribute 'UMAP'

From inside `spyder` the same script runs perfectly, but somehow when running it from the terminal it is not successful. I am new to python, so maybe I am not pointing correctly to `umap` library? Maybe it should be imported in a different way? Here is the whole script:

```
import umap
import numpy as np

def getopts(argv):
opts = {}  # Empty dictionary to store key-value pairs.
while argv:  # While there are arguments left to parse...
    if argv[0][0] == '-':  # Found a "-name value" pair.
        opts[argv[0]] = argv[1]  # Add key and value to the dictionary.
    argv = argv[1:]
return opts

if __name__ == '__main__':
  from sys import argv
  myargs = getopts(argv)
  print(myargs)
  path_to_save = myargs['-path']
  path_to_data = myargs['-data']
  data = np.loadtxt(open(path_to_data, "rb"), delimiter=",", skiprows=1)
  result = umap.UMAP().fit_transform(data.T)
  np.savetxt(path_to_save + "/umap.csv", result, delimiter=",")

```

Python version that is used by `spyder` and in the terminal seems to be the same: `3.6`. I am using `Mac OSX`.

# Answer
> 5 votes
You answered your own question, yet it is not a coincidence that you have to change the name of your script.

When Python executes your custom script, it reads:

> import umap

Python will first search the local sources (the directory of the script) , and afterwards the installed sources.

In this case Python identifies the local source to be `umap.py`, and imports the (currently executed) script, while foregoing the installed sources. Since you are calling functions, not defined in your current script, Python will not be able to fully execute your script.

For more information on importing see Module search path at docs.pyhton.org

# Answer
> 0 votes
It is pretty weird but the only thing that I needed to do is to rename my `python` script from `umap.py` to `umap_gen.py` and right away it started to work. So, I guess there can be such name collisions which we should avoid. Although I would prefer the file to be named as `umap.py` still...

---
Tags: python, umap
---


# Question
Title: Can you query the RCSB PDB API with a ligand and return the IDs of all entries in which it appears as a free ligand?
If you're looking at the entry for a ligand, for example, adenosine, there's a link on the webpage to entries for which this entry exists as a free ligand. However, I want to access this via the API and can't seem to find docs/info on how to do that.

I have been looking at the fields you can use to get custom reports.

# Answer
> 2 votes
I'm not sure about the RCSB, but you can get such a list from the PDBe REST API, as documented in the compounds section:

http://www.ebi.ac.uk/pdbe/api/pdb/compound/in\_pdb/ADN

Alternatively, from an older service:

http://www.ebi.ac.uk/pdbe-srv/pdbechem/PDBEntry/download/ADN

---
Tags: pdb
---


title: Violinplot of gene expression

text:
# Question
Title: Violinplot of gene expression
I have plotted the log normalized expression of two genes by violonplot for 4 clusters. I have links to my pictures and Seurat object too. My problem is this; in violin plot I can not see the mean or any centennial tendencies so that I don't know if two genes is expressing higher or lower in contrast to each other in each cluster.

violonplot vlnPlot(seuset,c(""))

Actually when I am thinking I see I need a violin or box plot showing the ratio of these genes to each other in each cluster and how nice to have mean also in the plot

Somebody knows how I can do that?

# Answer
You want to (1) see the mean for each gene, and also to (2) calculate a ratio of expression levels of two genes, then compare it between clusters.

(1) First, notice that `vlnPlot()` is deprecated. Use `VlnPlot()`.

You can try using the parameter `do.sort=T`:

`VlnPlot(object=seuset, features.plot=c("DDB_G0267412", "DDB_G0277853"), do.sort=T)`

Alternatively, you can return the ggplot2 object and then plot the means etc:

```
p <- VlnPlot(object=seuset, features.plot="DDB_G0277853", do.return=T)
p <- p + geom_boxplot(width=0.05)
p

```

---

(2) There are a few problems with calculating the ratio of gene expression levels. Firstly, what do you mean by gene expression level and how do you measure it? Surely you can have the read counts, but how do you interpret them? Is 30 counts a high level? Is it much more than 60 counts, or is it roughly the same? The same applies to the calculated ratios and the differences between them, even if we ignore amplification, gene length and other biases. And it is very hard to interpret ratios if the reference can also change. The problem is somewhat similar to RT-qPCR, where people use a set of reference genes whose expression has previously been shown to be invariant under the conditions.

In theory, you could use the raw counts (object@raw.data), the log + normalized counts (object@data), or the scaled counts (object@scale.data). The raw counts are biased by sequencing-depth, and the ratio of log or scaled values are not easily interpretable or intuitive. Therefore the library-size normalized (non-log) values seem to be the best.

Additionally, you could calculate the ratio of two genes either (a) for each cell (paired), or (b) for each group. (a) is problematic, because of the zero values: you will have many NaN and Inf values, which cannot be removed without biasing the data.

Here are some functions for retrieving and plotting data from the object:

```
ident1 <- WhichCells(seuset, ident=1)
ident2 <- WhichCells(seuset, ident=2)

FetchData(seuset, vars.all='DDB_G0267412', cells.use=ident1, use.raw=T)

datatoplot1 <- seuset@raw.data['DDB_G0267412', ident1] / seuset@raw.data['DDB_G0277853', ident1]
datatoplot2 <- seuset@raw.data['DDB_G0267412', ident2] / seuset@raw.data['DDB_G0277853', ident2]

boxplot(datatoplot1, datatoplot2)

stripchart(datatoplot1,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=1)

stripchart(datatoplot2,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=2)

```
> 4 votes

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


title: Spearman correlation for large dataset

text:
# Question
Title: Spearman correlation for large dataset
I have two datasets (DataA and DataB) and I want to find the Spearman correlation between genes and also pull out the gene names (stored in first column of dataset) in `R`. I am using `fread` from `read.table` to read the file and `cor.test` to find `Rho` and `p-value`. The code works fine for a smaller data set of 10 obs, 10 variables. But, my actual data set is larger (DataA - 795 obj, 542 variables and DataB - 925 obs, 542 variables). When I try to run the code on my actual data, it is consuming a large amount of time (it has been 23 hrs since my code is running). Is there a way to optimize this code? Or is there a bug in the code?

```
library(data.table)
A.data <- fread("DataA.csv", header = TRUE, data.table = FALSE, stringsAsFactors = FALSE)
B.data <- fread("DataB.csv", header = TRUE, data.table = FALSE, stringsAsFactors = FALSE)

correlation.result <- data.frame()

## Computing Correlation

for (i in 1:nrow(A.data)) {
  for (j in 1:nrow(B.data)) {
       correln <- cor.test(as.numeric(A.data[i,2:ncol(A.data)]), as.numeric(B.data[j,2:ncol(B.data)]), method = "spearman", exact = FALSE)
       rho <- correln$estimate
       p.val <- correln$p.value
       Gene <- paste(A.data[i,1], B.data[j,1], sep="_", collapse=" ") # To pullout gene name pairwise
       comb.data <- data.frame(rho, p.val, Gene)
       correlation.result <- rbind(correlation.result,comb.data)
 }
}

```

# Answer
> 3 votes
Your code looks fine, but `rbind` to progressively build and resize a data frame can be gluttonous in terms of memory and time, especially given the A-row x B-row times you are doing this operation, which is fine for very small datasets, but not often realistic for real-world data.

Instead, perhaps `cat` results to a text file, where results are written in tab-delimited columns. You can then read in that text file later on.

In addition to serializing or saving your results to text, which make them easier to work with in non-R contexts, you won't have to deal with R's memory allocation problems when running `rbind` repeatedly. You just write out the answer to a file and read it in later on, in one pass.

Another approach is to preallocate memory to a presized data frame, and then write to that allocated space.

You know the size of your inputs ahead of time, which makes this approach tenable. In this case, it looks like you are storing two floats and a string (`rho`, `p.val`, and `Gene`) and that you are going to store A-row x B-row number of them.

Here's one way you might go about this:

```
# create a presized data frame
Arows <- nrow(A.data)
Brows <- nrow(B.data)
correlation.result <- data.frame(t(matrix(c(0.0, 0.0, "foo"), nrow=3, ncol=Arows * Brows)), stringsAsFactors=FALSE)
colnames(correlation.result) <- c("rho", "p.val", "Gene")

# assign values to that presized data frame
Acols <- ncol(A.data)
Bcols <- ncol(B.data)
for (i in 1:Arows) {
  # get A data
  Ai <- as.numeric(A.data[i, 2:Acols])
  Aname <- A.data[i, 1]
  for (j in 1:Brows) {
    # calculate row offset in presized data frame
    ij <- (i - 1) * Brows + j
    # get B data
    Bj <- as.numeric(B.data[j, 2:Bcols])
    Bname <- B.data[j, 1]
    # run test
    correln <- cor.test(Ai, Bj, method = "spearman", exact = FALSE)
    rho <- correln$estimate
    p.val <- correln$p.value
    Gene <- paste(Aname, Bname, sep="_", collapse=" ")
    # assign value to row of result data frame
    correlation.result[ij,] <- c(rho, p.val, Gene)
  }
}

```

This avoids R needing to resize and recopy the data frame contents as they grow, which I suspect is the major reason it is taking so long to run.

Chapter or "Circle" 2 of Patrick Burns' book *The R Inferno* discusses the problems with growing objects with `rbind` and ways to handle it — you might find this a useful read. This book is a PDF download via your search engine of choice.

# Answer
> 1 votes
You can use `corr.test()` from psych package, e.g.

```
ct <- corr.test(t(A.data), t(B.data), method='spearman', adjust='none')

```

then `ct$p` and `ct$r` would give you (uncorrected) p-value and correlation matrices of size `795 x 925`. Note that `t()` is required to compare matrices row-wise.

That would be orders of magnitude faster than loop-based approach.

---
Tags: r
---


# Question
Title: Can I ask DESeq2 which variable alone explains which gene's behavior the best?
Imagine I have a dataset of highly periodic gene expression, taken at densely spaced timepoints. I know that the expression is periodic, but different genes have different offsets.

It is trivial to convert my timestamps into a cosine (`cos_time = cos(offset + time*pi*2/PERIOD)`), so that this new variable is linear with the genes' expressions -- but only with the ones that have a peak at `offset`.

What if I could have multiple `cos_time`s (let's call them `phases`) and then ask DESeq2: "OK, which gene's behavior is most linear with which phase?" Barring comparing p-values of each LRT test, is there a better way?

Disclaimer: *I know about cosinor, Lomb-Scargle periodograms, ARS, JTK etc, and I know they are likely more applicable in this particular case. I'm interested in whether it is in general possible to ask such a question from DESeq2, even if the variables in question are not periodic. Phase just seemed like the simplest concept to demonstrate this on.*

# Answer
> 2 votes
If you have densely spaced time points the using non-linear regression to directly fit the signal to a cosine or sine curve, which would directly produce the offset, would be much simpler. DESeq2 wouldn't be the appropriate tool then, you'd use base R and the `nls()` function.

# Answer
> 1 votes
If I were to do this, my first guess would be to add in the phases as a separate factor in the statistical model, something like:

```
~ Condition + sex + phase1 + phase2 + phase3 + ...

```

---
Tags: r, rna-seq, differential-expression, deseq2
---


title: Better aligner than bowtie2?

text:
# Question
Title: Better aligner than bowtie2?
`Bowtie2` is probably the most widely used aligner because of it's speed. Burrow-wheeler (BW) algorithms (including `bwa`) tend to be faster. However, they have limitations when it comes to aligning very short reads (e.g. gRNA). Also, setting maximum number of mismatches allowed is complicated by the seed length, overlaps and other parameters.

I wonder if there is any better multi-purpose aligner out there. May be with algorithm other that BW. One which allows special cases e.g. allowing shord reads and high number of mismatches.

# Answer
> 16 votes
Bowtie2 is no longer the fastest aligner. Salmon and Kallisto are much faster, but have been designed to optimise RNASeq mapping. Their speed is gained from avoiding a strict base-to-base alignment, but they can output mostly-aligned reads (i.e. position-only, without local alignment) as pseudo-alignments. See here for more details.

Both Kallisto and Salmon can do additional bootstrapping that is interpreted by *sleuth* for improved performance in isoform detection. They can also output counts that are equivalent to read-level counts from other programs, which can then be used by other downstream gene-based differential expression analysis software (e.g. DESeq2). Salmon has additional options that can correct mappings for sequence-level and GC bias.

HISAT2 is from the same group as Bowtie2, and does the same sort of stuff, but with a few optimisations added on top. In particular, it's much better at working out split reads from RNASeq runs. Like Bowtie2, it will do local alignment of reads.

For quick alignment of long reads, minimap2 works well. For high-accuracy alignment (but comparatively slower), LAST works well.

Most bioinformaticians seem to prefer STAR for things that Bowtie2 was previously used for. I'm not yet convinced it's a better alternative, and currently prefer HISAT2 for high accuracy short-read alignment.

---
Tags: sequence-alignment, sequence-analysis
---


# Question
Title: CNSeg no method coercing S4 class to vector
I'm attempting to analyze some genomic data with the help of a package called CNSeg. The package is located here and the documentation is here.

My input data looks like:

```
> str(segments)
'data.frame':   11897 obs. of  7 variables:
 $ X                : int  0 1 2 3 4 5 6 7 8 9 ...
 $ SegmentID        : int  72 73 74 75 76 77 78 79 80 81 ...
 $ Chromosome       : int  1 1 2 2 2 3 3 3 3 3 ...
 $ StartPosition    : int  754192 145260908 21494 141215321 141275624 63411 69812903 69884262 126473310 126790130 ...
 $ StopPosition     : int  145258178 249212878 141214996 141275051 243052331 69811900 69884106 126457276 126772699 197852564 ...
 $ Median.Log2.Ratio: num  -0.014 0.311 -0.003 0.059 -0.012 -0.018 -0.106 0.007 -0.171 0.001 ...
 $ FileName         : Factor w/ 95 levels "TSB02183","TSB02184",..: 1 1 1 1 1 1 1 1 1 1 ..

```

(Note: I did import the segments data with FileName as a vector instead of a factor as well. No difference in the error below).

I followed the vignette to the letter and, at the end of the vignette, the creator tells me that I can use the dist() method. When I try to use it, I get an error: **Error in as.vector(data) : no method for coercing this S4 class to a vector**.

Here is my session information:

```
> segments <- read.csv("Probe_Segments_CN.csv")
> cnseg <- CNSeg(segList = segments, chromosome = "Chromosome", end = "StopPosition", start = "StartPosition", segMean = "Median.Log2.Ratio", id = "FileName")
> rdseg <- getRS(cnseg, by = "region", imput = FALSE, XY = FALSE, what = "mean")
Processing samples ... Done
> data("geneInfo")
> geneInfo <- geneInfo[sample(1:nrow(geneInfo), 2000), ]
>
> rdByGene <- getRS(cnseg, by = "gene", imput = FALSE, XY = FALSE, geneMap = geneInfo, what = "median")
>
> reducedseg <- rs(rdseg)
> f1 <- kOverA(5, 1)
>
> ffun <- filterfun(f1)
>
> filteredrs <- genefilter(rdseg, ffun)
> filteredrs <- madFilter(rdseg, 0.8)
> dist(filteredrs)
Error in as.vector(data) :
  no method for coercing this S4 class to a vector
> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] tools     stats     graphics  grDevices
[5] utils     datasets  methods   base

other attached packages:
[1] CNTools_1.34.0    genefilter_1.60.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.17         AnnotationDbi_1.40.0
 [3] BiocGenerics_0.24.0  splines_3.4.0
 [5] IRanges_2.12.0       bit_1.1-14
 [7] lattice_0.20-35      xtable_1.8-2
 [9] blob_1.1.1           parallel_3.4.0
[11] grid_3.4.0           Biobase_2.38.0
[13] DBI_1.0.0            survival_2.41-3
[15] bit64_0.9-7          digest_0.6.15
[17] Matrix_1.2-9         S4Vectors_0.16.0
[19] bitops_1.0-6         RCurl_1.95-4.10
[21] memoise_1.1.0        RSQLite_2.1.1
[23] compiler_3.4.0       stats4_3.4.0
[25] XML_3.98-1.11        annotate_1.56.2

```

# Answer
I was not able to get the CNSeg clustering method to work, but the work-around was to export the data and cluster it independently of the CNTools package (thanks to Llopis's answer below). I'll post the analysis below.

---

```
require(CNTools)
segData <- read.csv("result_cnv.csv", stringsAsFactors = FALSE)
head(segData)

# Create inital CN object
cnseg <- CNSeg(segList = segData, chromosome = "Chromosome", end = "StopPosition", start = "StartPosition", segMean = "Median.Log2.Ratio", id = "FileName")
cnseg

# Create inital RD object
rdseg <- getRS(cnseg, by = "region", imput = FALSE, XY = FALSE, what = "mean")
rdseg

# Collect gene information
data("geneInfo")
geneInfo <- geneInfo[sample(1:nrow(geneInfo), 2000), ]

# Create an RD based on gene information
rdByGene <- getRS(cnseg, by = "gene", imput = FALSE, XY = FALSE, geneMap = geneInfo, what = "median")

# Initalize reduced segment
reducedseg <- rs(rdseg)

# Create a function that evaluates to TRUE if at least 5 of the argument elements are larger than 1
f1 <- kOverA(5, 1)

# Create a filter based on f1
ffun <- filterfun(f1)

# Use the CNTools genefilter
filteredrs <- genefilter(rdseg, ffun)

# Use the CNTools madFilter
filteredrs <- madFilter(rdseg, 0.8)

filteredrs

# Write the filtered data to file
write.csv(attributes(filteredrs)[1], "CNseg_CNOut.csv")

# Reading the recently-written file helps with formatting
CNseg <- read.csv("CNseg_CNOut.csv")

# Keep only the columns with actual sample values
CNseg <- CNseg[5:ncol(CNseg)]

# Calculate eluclidian distance between samples
d <- dist(t(CNseg), method = "euclidean")

# Calculate clustering
hc1 <- hclust(d, method = "ward.D")
plot(hc1, cex = 0.6, hang = -1, main = "Clusters of Copy Number Alterations", xlab = "Euclidean Distance")

```
> 1 votes

# Answer
If you are using the latest version of Bioconductor you should always look at the latest documentation of the package in Bioconductor. There could be changes from 2011 to 2018 (In this case there aren't any changes in the section you mention).

The problem is the package doesn't define how `dist` should behave with that object class (filteredrs is a RS object if I read the source code well).

You can extract the metrics to calculate the distance yourself or ask the maintainer to add this method to the class. However it could be that the maintainer answer you if you post the question in the support site of Bioconductor.
> 2 votes

---
Tags: r, bioconductor
---


# Question
Title: Why does `pratchet()` in `phangorn` not find the most parsimonious trees?
I called `pratchet()` on a dataset with binary values for the characters and it returned a parsimony value of 35. I then calculated the parsimony value for the dataset with a tree that I myself constructed and it was 7. How is this possible at all? I thought that `pratchet()` was supposed to return the most parsimonious tree on the basis of the data set, but that clearly did not happen in this case. Is there other software that is guaranteed to find the most parsimonious tree?

# Answer
Perhaps the tree which you have calculated is not being caught in the search heuristic? Have you tried (as a trial) setting the tree which you calculated as the starting topology, just to check that it is indeed being calculated on by pratchet()?

Alternatively, you could try playing around with the search heuristic settings, like varying from ratchet to NNI, or you could also try using an external program, like PAUP\*.
> 1 votes

# Answer
For dendrograms with less than about 10 OTUs/taxa/leaves it is possible to evaluate all possible dendrograms using an exhaustive search algorithm. However, for more than about 10 OTUs, heuristics need to be employed to try and "home-in" on and score a smaller set of "better" trees. However, these heuristics are not guaranteed to find ***the*** best dendrogram.

Without some specific code it is hard to do more than speculate as to the problem you are encountering, especially where you get an error using data you manually create.
> 1 votes

---
Tags: software-recommendation, phylogenetics, algorithms, phylogeny
---


title: What are RIKEN genes?

text:
# Question
Title: What are RIKEN genes?
I am using data from the Tabula Muris Consortium.

These are gene counts from scRNA-seq of mouse cells.

There are some specific genes with name ending with suffix 'Rik' (e.g. 0610005C13Rik or 0610007C21Rik). Searching online I found that they are from RIKEN project.

However, they are different in some respects (e.g. antisense lncRNA gene, or mouse 0610007C21RIK Adenovirus).

I want to exclude them from downstream analysis, but it is unclear to me the main purpose of this collection. For which kind of analysis are they used?

# Answer
My understanding is that these are cDNA sequences from the RIKEN project that don't yet have another gene name. The project carried out a comprehensive sequencing of the mouse genome, and as a result found a lot of cDNA transcripts that were previously unknown. Now they have some annotation (which can be varying degrees of functional annotation), and an accepted official identifier to allow people to talk about these transcripts when discovered in other studies.

It'd be interesting to see if any novel transcripts have been found from the Tabula Muris project. They would also need their own identifier, possibly leading to similar intrigue about their history.
> 3 votes

# Answer
To extend a bit more to @gringer answer. RIKEN is a large Japanese research institute, that participates in different scientific collaborations.

One of these is the FANTOM collaboration (Functional ANnotation Of the Mammalian genome), which lately have resulted in the release of FANTOM CAT, a browsing tool for exploring lncRNA relative to associated Genes, Traits and Ontologies.

Many of the transcripts and genes discovered during this work, are not yet named, and are therefore called CAT-genes, following numbering comparable to ensembl (eg. CATG00000118335.1).

I'm unsure how RIK genes have been discovered, as I have not personally encountered these genes yet, but in case you likewise encounter more genes discovered in collaboration with RIKEN, you now also know of CAT genes.
> 2 votes

---
Tags: scrnaseq
---


# Question
Title: Is it possible to identify cells that are expressing two or more genes in Seurat?
I'm interested in looking into cells that are positive for two (or in some cases more) genes. I know I have some double positive just by looking at the FeaturePlot of those genes, but now I'm trying to figure it out how many are double positive and further percentage wise. Example:

50 cells are positive for Myh11
100 cells are positive for Acta2
30 cells are double positive
60% of Myh11+ cells are Acta2+
30% of Acta2+ cells are Myh11+

Is the data used to generated the FeaturePlot the scale.data? If so, I could just make a new dataset out of that, but I'm not sure if that is the case.

# Answer
The counts stored in the Seurat object are: raw counts (seuratobject@raw.data), the log + normalized counts (seuratobject@data), and the scaled counts (seuratobject@scale.data). `FeaturePlot()` plots the log + normalized counts.

In order to identify double-positive cells, you need to identify cells that express a gene (i.e. positive for a gene) and that is not trivial: you need a cutoff value for each gene, and also to distinguish between true and dropout zero counts.

Provided you solve these, you can retrieve your cells by:

```
MYH11.cutoff <- 1
ACTA2.cutoff <- 1
length(which(FetchData(seuratobject, vars.all='MYH11') > MYH11.cutoff))
length(which(FetchData(seuratobject, vars.all='ACTA2') > ACTA2.cutoff))
length(which(FetchData(seuratobject, vars.all='MYH11') > MYH11.cutoff & FetchData(seuratobject, vars.all='ACTA2') > ACTA2.cutoff))

```
> 4 votes

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


# Question
Title: Is removing samples based on clustering for downstream analysis a right choice?
I'm using TCGA Lung cancer data. I'm interested in doing differential analysis between Lung vs Normal. Before DEA, to check the distance between each pair of samples I plotted an MDS plot:

In this, I see some Tumor samples are clustered with Normals. As this is a two dimensional plot, it is a little hard to say.

So, I made a clustering heatmap with the top 10% genes:

(Where brown color is Normal and black color is Tumor). I also see some Normals clustering with Tumor here, and some Tumor are with the Normals.

My question is: Should I exclude those tumor samples which are clustering with the normals for differential analysis? Is it a good idea to remove samples based on clustering?

Note: I didn't see any technical bias of those samples

# Answer
I would be very hesitant to blindly exclude those samples based on the clustering. Check to see if the clusters actually denote some sort of batch effect, since it's not like all of the TCGA datasets were processed at the same place or time. Check for clinical covariates too. Also, do a PCA and check the genes with high loading on the relevant PC. If it appears that the genes are indicators for cancerous state, then it's likely the samples are just mislabeled and you can reasonably exclude them. Otherwise keep them in so your results better represent the underlying biological variability.

BTW, your heatmap and MDS plots are what I'd expect from a cancer dataset. You occasionally see a tumour with the normal samples and vice versa, since the problematic genes aren't always in the top 10% and there is likely some heterogeneity in the cancer samples.
> 5 votes

---
Tags: rna-seq, differential-expression, clustering
---


# Question
Title: determine if ChIP-seq peaks are broad or narrow
Is there a method to determine if the peaks are broad or narrow? ENCODE provides some guidelines:

Although those cover common histone marks, there are many others. If you are using one of the ones that are not listed, is there a way to determine the appropriate peak width?

This question was inspired by this previous post.

# Answer
I'd say mostly this is a question of understanding the underlying biology and the relevant literature.

If it is not known in the literature whether a mark is peaky or broad, evidence might come from FISH studies or Low throughput qPCR. Another way to look would be to examine the signal expressed as fold enrichment over input on a genome browser and look if peaks tended to be small (a few hundred bp, e.g. a transcription factor), or large (several 1kb, e.g. H3K27me3) or contain both (e.g. RNA polymerase II).
> 3 votes

---
Tags: chip-seq, peak-calling
---


# Question
Title: Where can I find the DNA sequence of the C. elegans rRNA cluster?
Where can I find the DNA sequence of the *C. elegans* rRNA cluster?

Ideally, I would be able to do it using the command-line.

I've tried using NCBI but the results are unsatisfactory. I am interested in the DNA sequence that is in the current reference genome of WormBase.

I have tried using the Bioconductor biomaRt package and create a table containing the description but when I grep 'rRNA', I only get proteins doing things with rRNA. I used the following code:

```
ensembl <- biomaRt::useMart('ensembl', dataset = 'celegans_gene_ensembl')
wormGenes <- biomaRt::getBM(attributes = c('external_gene_name', 'chromosome_name', 'start_position', 'end_position', 'gene_biotype', 'description'), mart = ensembl)
wormGenes[grep('rRNA', wormGenes$description), 6]

[1] "rRNA 2'-O-methyltransferase fibrillarin  [Source:UniProtKB/Swiss-Prot;Acc:Q22053]"
[2] "Pre-rRNA-processing protein pro-1  [Source:UniProtKB/Swiss-Prot;Acc:Q22006]"
[3] "Pre-rRNA-processing protein TSR1 homolog  [Source:UniProtKB/Swiss-Prot;Acc:Q19329]"
[4] "rRNA methyltransferase 2, mitochondrial  [Source:UniProtKB/Swiss-Prot;Acc:O62251]"
[5] "Probable rRNA-processing protein EBP2 homolog  [Source:UniProtKB/Swiss-Prot;Acc:Q09958]"
[6] "Putative rRNA methyltransferase  [Source:UniProtKB/TrEMBL;Acc:Q5CCJ2]"

```

# Answer
The info is actually in the `gene_biotype` column of my `wormGenes` data frame. I could select the rRNA genes using `dplyr::filter(wormGenes, gene_biotype == 'rRNA')`. Unfortunately, no description is available for the genes, so I will have to get it using the gene IDs.
> 1 votes

---
Tags: database, biomart, data-download, sequence-annotation
---


title: Input normalization in ChIP-seq

text:
# Question
Title: Input normalization in ChIP-seq
If I subtract input counts from ChIP counts (for every gene, since I have one peak per gene) I get negative values for most genes. This makes sense to me, because (as can be seen in the figure) input and IP have same sequencing depth but IP is very biased towards few genes with lots of counts.

Elsevier was informed that the business model which they offered does not meet the demands. How should I implement input normalization?

P.S. The subtraction idea is from Bioconductor

# Answer
> 4 votes
I strongly suggest that you not try to come up with your own package for this when things like `CSAW` already exist in bioconductor and provide a number of useful normalization options.

For visualization purposes, I think it's best to simply take the log2 ratio (e.g., using `bamCompare` from deepTools) or alternatively normalize by total coverage (possibly excluding blacklist regions, if applicable for your organism).

# Answer
> 1 votes
If you want to compare conditions (or strains), I don't see why you need Input at all (unless you want to define intervals).

It isn't clear from you question (and subsequent comment) on what you want to achieve, but I can see 2 possible scenarios:

1. **You already have predefined intervals to test on**

In your case (Pol2-ChIP), this might be for example gene promoters. You could compare the Pol2-ChiP libraries directly using standard test for count data as established for RNA-Seq (i.e. edgeR, baySeq etc) in which the intervals are genes (as opposed to promoters or so). No need for input in this scenario, all you need are ChIP libraries for both conditions (don't forget to do replicates).
2. **You don't have predefined intervals**

You could defined them through **peak calling**, in this case you'd in fact need input. A nice approach for this is described in Diaz et al and implemented in deeptools. It poses an elegant method to separate the IP signal from background input, which is probably what you want I suppose. Having defined peaks, you could do scenario 1) using them as intervals.

---
Tags: bioconductor, normalization, chip-seq
---


title: KEGG FTP vs KEGG API

text:
# Question
Title: KEGG FTP vs KEGG API
I was reading the KEGG plea and I found that it doesn't forbid using the KEGG API. Then, what is in the FTP server license for personal use/academic use that it is not covered by the API?

Or I could download all the database via the API?

PS: I couldn't find how to download which genes are in which pathways but it might be I didn't expend enough time with the queries.

# Answer
I found the way to download the genes and pathways thanks to this answer using KEGGREST.

It seems that the API allows to fully download the database because it allows to download info, list, find, get, conv, link, ddi.

In fact in the message of 2011-05-21 it mentions that the API includes more options to download the KGML files that previously where only possible to get via the FTP site.
> 2 votes

# Answer
If you like to get pathways and corresponding genes as simple table format in a plain text file, use GeneSCF 'prepare\_database' module. This can be simple solution.

**FTP vs API**:

The differences I see is,

* quick access to the updated content via FTP. With that I assume using FTP version you will have access to weekly database updates. Also bulk download is possible and more convenient via FTP subscription.
* But KEGG REST API updates once in month or two months or a year (not a recently updated contents).
> 1 votes

# Answer
We implemented Bio2BEL KEGG package to download the KEGG database via its RESTful API
> -1 votes

---
Tags: database, kegg
---


# Question
Title: Obtaining uniquely mapped reads from BWA mem alignment
*This is based on a question from betsy.s.collins on BioStars. The original post can be found here.*

> Does anyone have any suggestions for other tags or filtering steps on BWA-generated BAM files that can be used so reads only map to one location? One example application would be to find seeds for the TULIP assembler/scaffolder, which works best for reads that map to unique genomic locations.

The questioner is referring to genomically-unique alignments, which are something different from multiple possible alignments at a single location. If the only multiple alignment is at the same genomic location, those alignments are still unique in the genomic sense. The impossibility of finding the exact, correct, alignment is a well-known problem, and there are a few downstream methods (e.g. left normalisation) to make sure that multiple local alignments and/or sequencing errors have reduced impact.

The concept of "uniquely mapped reads" is a loaded term, and most sources suggest filtering by MAPQ should do the trick. However, this approach doesn't seem to work when using BWA as a read mapper.

Uniquely mapped reads (i.e. a read that maps to a single location in the genome) are sometimes preferred when running analyses that depend on the quantification of reads, rather than just the coverage (e.g. RNASeq). Duplicated reads require additional filtering or normalisation on the analysis side, and most downstream programs won't consider concepts like "a fifth of the read maps somewhere here on chromosome 5, and four-fifths of the read maps somewhere here on chromosome 14."

# Answer
**Update** \- as of January 2021, samtools can now do filtering based on an expression that includes tag variables. In this case, this expression can be used to exclude any reads that have either an XA or SA tag:

```
samtools view -b mapped.bam -e '!([XA] | [SA])' > unique_mapped.bam

```

For more details on the samtools expression parameter, see the samtools documentation:

https://www.htslib.org/doc/samtools.html#FILTER\_EXPRESSIONS

---

Original answer follows....

To exclude *all possible* multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use `grep` on the uncompressed SAM fields:

```
samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique_mapped.bam

```

Explanation follows...

I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM file produced by BWA, and has no way of getting the original reads. One high-effort solution would be to extract the mapped reads from the BAM file and re-map with a different mapper that uses the MAPQ score to indicate multiple mappings.

... but what if that were not possible?

My understanding of BWA's output is that if a read maps perfectly to multiple genomic locations, it will be given a high mapping quality (MAPQ) score for both locations. Many people expect that a read that maps to at least two locations can have (at best) a 50% probability of mapping to one of those locations (i.e. MAPQ = 3). Because BWA doesn't do this, it makes it difficult to filter out multiply-mapped reads from BWA results using the MAPQ filter that works for other aligners; this is likely to be why the current answer on Biostars \[`samtools view -bq 1`\] probably won't work.

Here is an example line from a BWA mem alignment that I've just made. These are Illumina reads mapped to a parasite genome that has a *lot* of repetitive sequence:

```
ERR063640.7     16      tig00019544     79974   21      21M2I56M1I20M   *       0       0       TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13      AS:i:77 XS:i:67 XA:Z:tig00019544,-78808,21M2I56M1I20M,6;tig00019544,-84624,79M1I20M,6;tig00019544,-79312,33M4I42M1I20M,8;

```

BWA mem has found that this particular read, ERR063460.7, maps to at least three different locations: tig00019544, tig00019544, and tig00019544. Note that the MAPQ for this read is 21, so even though the read maps to multiple locations, MAPQ can't be used to determine that.

However, the alternative locations are shown by the presence of the `XA` tag in the custom fields section of the SAM output. Perhaps just filtering on lines that contain the `XA` tag will be able to exclude multiply-mapped reads. The `samtools view` man page suggests that `-x` will filter out a particular tag:

```
$ samtools view -x XA output.bam | grep '^ERR063640\.7[[:space:]]'
ERR063640.7 16  tig00019544 79974   21  21M2I56M1I20M   *   0   0   TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13  AS:i:77 XS:i:67

```

... so it filtered out the tag (i.e. the tag no longer exists in the SAM output), but not the *read*. There are no useful bits in the FLAG field to indicate multiple genomic mappings (which I know can be filtered to exclude the read as well), so I have to resort to other measures.

In this particular case, I can use `grep -v` on the uncompressed SAM output to exclude alignment lines that have the `XA` tag (and re-compress to BAM afterwards, just to be tidy):

```
$ samtools view -h output.bam | grep -v 'XA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | grep '^ERR063640\.7[[:space:]]'
<no output>

```

Hurray! reads filtered. As a little aside, this `grep` search has a fairly substantial computational load: it's looking for some string with the text `XA:Z:` somewhere in the line, and doesn't actually capture every situation. Some masochistic person might come along at a later date and decide that they're going to call all their reads `HAXXA:Z:AWESOME!:<readNumber>`, in which case a tweak to this grep search would be needed to make sure that there's a space (or more specifically a tab character) prior to the `XA:Z:`.

Now I do a check for *any* duplicated read names, just to be sure:

```
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
ERR063640.1194
ERR063640.1429
ERR063640.1761
ERR063640.2336
ERR063640.2825
ERR063640.3458
ERR063640.4421
ERR063640.4474
ERR063640.4888
ERR063640.49
ERR063640.4974
ERR063640.5070
ERR063640.5130
ERR063640.5300
ERR063640.5868
ERR063640.6116
ERR063640.6198
ERR063640.6468
ERR063640.6717
ERR063640.6797
ERR063640.7322
ERR063640.750
ERR063640.7570
ERR063640.7900
ERR063640.8115
ERR063640.8405
ERR063640.911
ERR063640.9206
ERR063640.9765
ERR063640.9986

```

Oh... damn. I wonder what they are:

```
$ samtools view output_filtered.bam | grep '^ERR063640.3458[[:space:]]'
ERR063640.3458  16  tig00002961 5402    60  58S38M  *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA    DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJDHLNKCK;KMKGMFKJILIIIMKI9JLKKHEJFII?CC    NM:i:0  MD:Z:38 AS:i:38 XS:i:0  SA:Z:tig00002377,202353,-,14M3I5M1I35M38S,19,5;
ERR063640.3458  2064    tig00002377 202353  19  14M3I5M1I35M38H *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA  DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJ  NM:i:5  MD:Z:5G48   AS:i:35 XS:i:27 SA:Z:tig00002961,5402,-,58S38M,60,0;

```

Aha! Supplemental alignments, which use the *official* `SA` tag \[other canonical alignments in a chimeric alignment\]. These appear to be situations where a single read has been split up, and maps to multiple locations. Note that in this case, ***both of the alignments still have MAPQ scores of over 3***. It sounds like the questioner would also want to get rid of these situations as well. This time, there are standard flag fields as well to deal with these situations (0x800: secondary alignment). Except it's not enough to just filter the supplemental alignment, because *both* read mappings should be removed, rather than just the one (or ones) that happened to be tagged as secondary.

Luckily, BWA appears to put the `SA` tag into all reads containing supplementary alignments (if this is not the case, I'm sure someone will correct me on that). So, I add in the `SA` search as an additional `grep` filter:

```
$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
<no output>

```

Done. Easy peasy! `</sarcasm>`

... that original "high-effort" solution of using a different aligner doesn't look so bad now.
> 13 votes

# Answer
For much faster and stable<sup>$</sup> implementation, sambamba can handle this using following one-liner:

```
sambamba view -t 12 -h -f bam -F "mapping_quality >= 1 and not (unmapped or secondary_alignment) and not ([XA] != null or [SA] != null)" test.bwa-mem.bam -o test-uniq.bam

```

See manpage and sambamba synatx for filters for details on further usage.

$: Besides `grep` based search being slow, it may return false ~~positive~~ negative hits (@gringer's reply above). I also do not find `awk` based approach reliable because fields like XA, SA, etc. are optional fields in SAM format, and not *positionally* fixed fields.
> 4 votes

---
Tags: sam, sequence-alignment, bwa
---


title: Public access to genomics databases

text:
# Question
Title: Public access to genomics databases
I'm a statistician, and am interested in applying the theory of algebraic Markov models to genomics. Here's one paper I'm interested in: Large-scale genome-wide association analysis of bipolar disorder identifies a new susceptibility locus near ODZ4. I'm also interested in other diseases (i.e. not just bipolar).

Are there any genomics databases that I can access publicly?

There are some links in Human Connectome Project, but I'm not sure if these are genomics databases (or which papers have used their research).

# Answer
> 3 votes
Most genomics data that is published should be available in GEO (or SRA). Usually in the paper they mention the accession number. However, not every study ends up in there, and sometime studies that are in GEO are not reproducible (poor description of the samples, or design). But if you look at the original papers, most of the times you'll find an accession number.

# Answer
> 0 votes
Another places you could look for SNP data would be 1000 genomes, the ExACT database and the genomAD database. Cancer specific data is available at the GDC and ICGC data portal.

Unfortunately a lot of what you'd like to work with is probably covered by restricted access due to privacy concerns. Lots of this data is stored in the dbGap database and the EGA. If you are an academic or employed by a company with a legitimate interest in the data you can usually apply for access to the restricted data, and I've never known anyone to be turned down if you have a legitimate reaons to access the data and can demonstrate your computer systems meet the required security standards.

# Answer
> 0 votes
Yet another place I find useful is provided by Repositive. You can think of it as a meta-search engine which includes most of the databases already stated and links hits to relevant papers.

With respect to your question about bipolar, this search could get you started.

---
Tags: database, human-genome
---


# Question
Title: What are some good sources of 3rd order gene expression data
I need good sources of gene expression data in the form of 3rd order tensors. Typically the commonly available datasets are in the form of a matrix, for instance, $sample \times gene$ or $gene \times time$. I need good sources for 3-dimensional data, for example $gene \times sample \times tissue$ or $gene \times sample \times time$ ... etc.

# Answer
I familiar with the inflammatory bowel disease which is quite complicated: It affects multiple sites and depending on when it appears it is different. So you can look at the GEO for microarrays and RNA-seq data of this disease and I'm sure you'll find datasets from the same patients from several tissues or time.

In this disease there is a multi'omics project that fits your needs. You can find (and download) data of high dimensionality: there are samples for the same time, tissue (but different locations), gene, over a period of time and for several omics (not only expression data).
The only problem is that there are some mistakes on the metadata.

However, I don't think they are usually described as 3rd order tensors and I doubt you'll find enough data to train a machine learning method.
> 1 votes

---
Tags: rna-seq, public-databases, microarray
---


# Question
Title: Human Genome Alignment and Variant Calling Benchmarks
I am doing a project that requires bench-marking the performance of computational tools for alignment and variant calling of human genome sequencing data. In particular I would like to establish benchmarks comparing the performance of GATK and ADAM/Avocado.

My main question is which genomic data-sets should I use as the input to the data processing pipeline. In Frank Nothaft's (author of ADAM and Avodaco) thesis he reports bench-marking the processing of NA12878 from the 1,000 genome's project in Avocado as compared to GATK's HalotypeCaller. He also reports benchmarking on 270 samples from the Simons Genome Diversity Data-set.

I am interested in reproducing this portion of his work, but do not know where I can get a hold of these datasets. Where can I get a copy of the NA12878 genome (BAM file) and also were can I get aligned samples from Simons Genome Diversity Data-set?

# Answer
> 2 votes
### NA12878

* My recommendation: download raw reads from Illumina BaseSpace and do alignment yourself. Aligning ~40X worth of human reads takes overnight with ~16 CPU cores.
* You can acquire Platinum Genome data from ENA. There is a BAM available for download, but that was done by Illumina's own aligner against a wrong reference genome. You'd better do alignment by yourself. Another problem with this dataset is it is too good to be representative. Modern data are worse.
* GIAB provides BAMs here:

```
ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/NIST_NA12878_HG001_HiSeq_300x/

```

It is not "typical", either.

### SGDP

* Instructions on downloading SGDP BAMs are available on this page, in particular this Word document. Warning: you need ~30TB disk space for these BAMs. I also heard it is tricky to download data.
* There is apparently a mirror from cancer genomics cloud. I don't know how easy it is to access.

---
Tags: gatk, adam
---


# Question
Title: Which data is being used for violin plot?
Sorry I just got totally confused conceptually.

If this is my raw count data

Likely, Seurat divides each value by sum of the column afterward times by 10000.

which gives so

```
                   s1.1        s1.2 s1.3       s1.4 s1.5
DDB_G0267178 0.00000000 0.009263254    0 0.01286397    0
DDB_G0267180 0.00000000 0.000000000    0 0.00000000    0
DDB_G0267182 0.00000000 0.000000000    0 0.03810585    0
DDB_G0267184 0.00000000 0.000000000    0 0.00000000    0
DDB_G0267188 0.02640801 0.000000000    0 0.01286397    0

```

Literally Seurat claims that final step would be taking natural log of the above matrix as I tested with

```
library(seurat)
mat <- matrix(data = rbinom(n = 25, size = 5, prob = 0.2), nrow = 5)
mat

    [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    1    1    1
[2,]    1    0    2    2    3
[3,]    1    1    1    1    2
[4,]    1    1    0    1    0
[5,]    0    0    0    1    3
mat_norm <- LogNormalize(data = mat)
mat_norm
[1,] 7.824446 .        7.824446 7.419181 7.014015
[2,] 7.824446 .        8.517393 8.112028 8.112028
[3,] 7.824446 8.517393 7.824446 7.419181 7.706713
[4,] 7.824446 8.517393 .        7.419181 .
[5,] .        .        .        7.419181 8.112028

```

My confusion is: seurat@data dose not give me natural log transformed data rather returns only data divided by column sum and scalled by 10000. So when plotting a violin plot if data are log normalised plot should not look as below picture

When I used

```
VlnPlot(object = y, features.plot = "DDB_G0277853", x.lab.rot = FALSE, y.log = TRUE)

```

plot changed totally. Could someone please help me in getting idea about seurat@data that is not really log transformed and by which data violin plot is being produced??? So if a gene has 10000 read counts, without log transformation violin plot should be so

This is my seurat object

# Answer
The normalized and log-transformed values are used for the violin plot. The argument `y.log` changes only the display of the data (scaling of the y axis).

Seurat has very good documentation. Section 7 in the FAQ explains what data is stored in the object:

> 7. How is data stored within the Seurat object? What is the difference between raw.data, data, and scale.data?
>
> **raw.data** The raw data slot (object@raw.data) represents the original expression matrix, input when creating the Seurat object, and prior to any preprocessing by Seurat.
>
> **data** The data slot (object@data) stores normalized and log-transformed single cell expression. This maintains the relative abundance levels of all genes, and contains only zeros or positive values. See ?NormalizeData for more information. This data is used for visualizations, such as violin and feature plots, most differential expression tests, finding high-variance genes, and as input to ScaleData (see below).
>
> **scale.data** The scale.data slot (object@scale.data) represents a cell’s relative expression of each gene, in comparison to all other cells. Therefore this matrix contains both positive and negative values. See ?ScaleData for more information.
> 2 votes

# Answer
Those two graphs don't look radically different. The y axis is definitely not log scaled on the top one, and it looks like it is scaled in the second one (y.log = TRUE)
> 2 votes

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


# Question
Title: How to find correlation between two specific genes in same dataset?
I would like to plot the correlation between two specific genes in my data. I have a matrix with genes in rows and samples in columns, with read counts data. If I want to look at the correlation between 2 genes

1) do I need to convert the counts to some other units?

2) Is normalization needed?

I currently have counts returned by `featureCounts`:

```
          Sample1  Sample2  Sample3  Sample4  Sample5

ESAM         803    338       528       841      712
ESCO1        1594   1401      847       392      611
ESCO2         274   778       288       77       204
ESD          6192   3502      2013      1573     1808
ESF1         1356   1497      1011      372      701
ESM1         398    58        114       582      135

```

I tried this way:

I transposed the above data making genes as columns and sample names as rows. I removed the rows. Lets say the genes as columns is in dataframe "df".

```
df.mat <- cor(df)
colnames(df.mat) <- NULL
rownames(df.mat) <- NULL
heatmap.2(df.mat, trace="none")

```

I have an error:

```
Error in .External.graphics(C_layout, num.rows, num.cols, mat, as.integer(num.figures),  :
  invalid graphics state

```

I want to look at the correlation between ESCO1 and ESCO2. Can anyone tell me how to do that and plot with a small example.

# Answer
> 3 votes
1. I suppose in an idealized world you'd use something like `salmon` to get TPMs, which can account for differences in effective gene length by sample. However that's quite likely to be overkill.
2. Not necessary, you're using a with-sample ratio, which wouldn't be changed by normalization.

Here is a small example in R, assuming your matrix or data frame is named `m`:

```
first = which(row.names(m) == "ESCO1")
second = which(row.names(m) == "ESCO2")
plot(m[first,], m[second,])

```

or `cor(m[first, ], m[second,])`. You'll want a rank-based correlation.

# Answer
> 0 votes
I'm answering my own question. I tried this way.

```
df <- read.csv("df.csv")

Samples ESAM    ESCO1   ESCO2   ESD ESF1    ESM1
Sample1 803     1594    274    6192 1356    398
Sample2 338     1401    778    3502 1497    58
Sample3 528      847    288    2013 1011    114
Sample4 841      392    77     1573  372    582
Sample5 712      611    204    1808  701    135

df2 <- data.frame(df[,-1],row.names = df[,1])

cormat <- round(cor(df2,use="pairwise.complete.obs"),2)
head(cormat)

library(ggcorrplot)
pdf("eg.pdf")
ggcorrplot(cormat,hc.order = TRUE,
           outline.col = "white",
           ggtheme = ggplot2::theme_gray,
           colors = c("#6D9EC1", "white", "#E46726"))
dev.off()

```

# Answer
> 0 votes
You can represent, and calculate, the correlations using psych:

```
library(psych)
pairs.panels(df,df,
             method = "pearson", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots
             ellipses = TRUE,xlim=c(0,500),ylim=c(0,500))

```

---
Tags: rna-seq, r, gene, bioconductor
---


# Question
Title: Correct way to assemble reference transcriptome - what's --samples_file option in Trinity?
I have 101 samples from 9 tissues of a nonmodel species, coming either from SRA or sequenced by my lab. I want to generate a reference transcriptome with Trinity for further analysis. In order to do this, I downsampled the reads in order to have 20 to 30M reads for each tissue, for a total of 270M PE reads, which I combined in two files, one for left reads and one for right ones.

Now, I came across the --samples\_file flag in Trinity and I'm not sure if I should use it instead of concatenating the transcriptomes. Any help on good practice on this?

# Answer
`--samples_file` flag is for DE analysis, you can put info about biological replicates in that file. You say that you want to make a ref transcriptome, so you don't need a samples file for DE analysis.

Why do you down sample? Are you not afraid you will not have enough reads to have a complete transcriptome now?
> 2 votes

---
Tags: assembly, transcriptome
---


# Question
Title: How to superimpose barplots in karyotype graphic using the chromPlot Bioconductor package?
I have a data frame containing 4 columns:

* Chrom: chromosome name of the feature
* Start: starting base of the feature
* End: end base of the feature
* Name: name of the feature, e.g. 5S ribosomal RNA

The head of the data frame is shown below:

```
      Chrom    Start      End                        Name
    1     1 15390242 15390720 16S ribosomal RNA (partial)
    2     1 15390900 15391794 16S ribosomal RNA (partial)
    3     1 15392798 15396046           23S ribosomal RNA
    4     1 15397439 15397917 16S ribosomal RNA (partial)
    5     1 15398097 15398991 16S ribosomal RNA (partial)
    6     1 15399995 15403243           23S ribosomal RNA

```

I want to plot these features on the karyotype of *C. elegans*. I am using the chromPlot package from Bioconductor.

I am using the following code:

* Subsetting the original data frame to get a data frame for each feature type: 5S, 16S and 23S rRNAs

```
annot1 <- dplyr::filter(cel, Name == '5S ribosomal RNA')
annot2 <- dplyr::filter(cel, Name == '23S ribosomal RNA')
annot3 <- dplyr::filter(cel, Name == '16S ribosomal RNA (partial)')

```

and then plot:

chromPlot(gaps = cel, figCols = 2, annot1 = annot1, annot2 = annot2)

This is OK because the two features are far away in the genome. However, if I try:

```
chromPlot(gaps = cel_gap, figCols = 2, annot1 = annot1, annot2 = annot2,
          annot3 = annot3)

```

one of the features masks the other because they are close to each other in the genome. How can I make the green bar stack on top of the yellow one instead of covering it?

# Answer
> 2 votes
```
chromPlot(gaps = cel_gap, figCols = 2, annot1 = annot1, annot2 = annot3, annot3 = annot2,
          chrSide=c(1, 1, -1, 1, -1, 1, -1, 1))

```

The `chrSide` argument decides on which side the histogram goes, so it can be used to solve my problem. The bars are not stacked but the result is the same: more accurate visualisation.

Thanks to Llopis for the guidance.

---
Tags: r, bioconductor, visualization, chromplot
---