# Question
Title: Is there a tool that can perform a read-group-aware mpileup from a single file?
I would like to perform a `samtools mpileup` from a single file that contains thousands of read groups with different SM tags. I could split the bam by read group using `samtools split`, and then perform the mpileup, but splitting into thousands of files and then performing an mpileup from thousands of files is crazy slow on my machine.

I could imagine an implementation using pysam, but I was hoping someone has a solution already.

Bcftools mpileup already has this feature, but the feature was never implemented for the original mpileup format, which is what I need.

### An example

Contents of `input.sam`:

```
@HD VN:1.6  SO:unknown
@SQ SN:chr1 LN:1000000
@RG ID:0    SM:sample_0
@RG ID:1    SM:sample_1
r0  0   chr1    24  0   1M  *   0   0   G   I   RG:Z:0
r1  0   chr1    24  0   1M  *   0   0   G   I   RG:Z:1

```

When I run samtools mpileup I get this behavior:

```
> samtools mpileup input.sam
[mpileup] 2 samples in 1 input files
chr1    24  N   2   ^!G$^!G$    II

```

Whereas, what I actually want is this output:

```
[mpileup] 2 samples in 1 input files
chr1    24  N   1   ^!G$    I   1   ^!G$    I

```

For a long time I thought this was a bug in samtools mpileup, but it appears that it is a feature that was described/hinted at in the samtools documentation, but was never implemented: https://github.com/samtools/samtools/issues/599#issuecomment-604941884

# Answer
> 2 votes
I broke down and finally scratched this itch. I have implemented a tool that performs read-group-aware text-pileup from a single SAM/BAM file. The tool is called streaming\_pileupy and it's available from pypi for installation with pip. After installation the command would be:

```
spileup input.bam sample_names.txt

```

The tool is pre-alpha and missing some features. Contributions are welcome.

---
Tags: sam, samtools, pysam, mpileup
---


# Question
Title: How to combine two Genome-wide Association Study (GWAS)?
I did a GWAS analysis in the past for antibiotic resistance of E. Coli and the results were interesting (matching the literature). I did a new GWAS analysis for some new samples, but the results are not in agreement with the previous GWAS and the literature. So, I was wondering if the combination of them would be in agreement with the literature, as having more samples to support an argument is always a good idea.

So, my question is how I can combine two GWAS study and what quality control steps do I need to follow?

# Answer
The first way is to combine two samples and re-run GWAS. A second way is to do a meta-analysis using GWAS results.
> 0 votes

# Answer
I'd recommend first having a think about how you would go about combining two GWAS with the *same* sample set (e.g. split the samples randomly in two). Even with the same source population, you'll get different associations coming up by sampling a different sub-set of that population.

Once you can do that multiple times, you can attach a level of confidence to each SNP association within each study, which can be treated as a prior probability for the purpose of combining multiple independent studies.

It sounds like you're in a great position to be able to do this, because you've done the analysis on both groups and \[presumably\] have all the genotype-level data from both studies.
> 0 votes

---
Tags: sequence-alignment, ngs, sequence-analysis, genomics, gwas
---


# Question
Title: What is outliers and how to determine outliers in the following phylogenetic tree?
I an confused with the term 'outlier' and 'outgroup' in phylogeny. Could someone please explain with the example phylogenetic tree above which is the correct term? Thank you very much !!!

# Answer
Distance- and ML-based algorithms using reversible models can't find the root of trees. A classical method to root a tree is to use an outgroup (not outlier), which is a species/sequence *known* to directly descend from the root of the tree. In your case, it is relatively easy: add a chicken or fish ortholog to your dataset and put the root on the chicken/fish branch.

Outgroup works well for single-copy genes used in species tree construction (e.g. ADH). However, it doesn't always work. The culprit is the "known" part. Say a gene has two copies A and B in vertebrates. Copy A was lost in rodent and copy B lost in primate. If you choose a chicken A gene as the outgroup, the correct tree should be ((primary-A,chicken-A),rodent-B). Without knowing the true history, you may forcefully put the root at the chicken-A branch and build a wrong tree ((primate-A,rodent-B),chicken-A).

There are a few other tree rooting methods. An easy approach is to put the root at the longest branch in the tree, assuming the presence of molecular clock. When the species tree is known, you can root a gene tree by minimizing the number gene duplication/loss events in the history. I generally prefer the latter approach when the relevant information is available.
> 3 votes

# Answer
you may use the OD-seq algorithm to filter outliers in a multiple alignment. https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0702-1
> 0 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: How do I pull singe cell RNA sequencing data from GEO database?
I am new to R and computational biology. I am trying to look through a published data set to check for gene expression for my own project. I am having trouble finding materials to teach me how to analyze single cell RNA sequencing data on R from the GEO NCBI database.

This is the data that I am trying to access https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138651

I would appreciate any directions or suggestions for tutorials.

Thanks!

# Answer
Install package

```
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GEOquery")

```

Use function getGEO to download the files.

```
expressionSet <- getGEO("GSE138651")
expressionSet <- expressionSet[[1]]

```

You can then look at the info of the samples:

```
pData(expressionSet)

```

And the assay data itself

```
exprs(expressionSet)

```
> -1 votes

# Answer
The simplest would be using a count matrix (at the end of the link you have shared, section "Supplementary file"). For example `GSE138651_Vagal_WholeNodoseSeq_raw.mtx.gz` corresponds to "raw" count matrices obtained after mapping and counting the reads. This file, along with the `GSE138651_barcodes.tsv.gz` and `GSE138651_genes.tsv.gz` can be used to "read" the data into `R`; to that end there are many options in the form of packages, Seurat is one of the best documented and hence is beginner-friendly.
> 5 votes

# Answer
You might want to check out the BioConductor *`GEOquery`* package.

Install using

```
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GEOquery")

```

More information on how to use this package and read in data from GEO can be found in the vignette
> -1 votes

---
Tags: r, scrnaseq, phylogenetics, geoquery
---


# Question
Title: How to remove sequence reads contains more than 2 X from multifasta file?
I have 5000 protein sequences in one multifasta file. I found more reads have gaps as X in their reads. So, want to eliminate those reads completely (Whole protein seq) from the file. I am keeping filter criteria as if a read contains morethan 2 X ( continuesly or anywhere in the read) should be removed. Thanks in advance for your help.

The input sequence looks like this

```
>Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF
>Prot2
ANSTVKKKKLLLYYYSSSEERXXXXXXXXXXXFGHYFGHYFGHFYVHFGFYVHCEDYHF
>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

```

I want output Like this

```
>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

```

# Answer
A Biopytnon solution that works for `.gz` files too :-). The script make a script `filter_faa.py`:

```
#!/usr/bin/env python3

import sys
from Bio import SeqIO
from mimetypes import guess_type
import gzip

filename = sys.argv[1]

encoding = guess_type(filename)[1]
_open = partial(gzip.open, mode='rt') if encoding == 'gzip' else open

with _open(filename) as ffile:
    for seq_record in SeqIO.parse(ffile, 'fasta'):
        if seq_record.seq.count('X') < 3:
            print(seq_record.format("fasta"))

```

Then usage is

```
python3 filter_faa.py my.faa > my_filtered.faa

```

The `my_filtered.faa` then will contain

```
>Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF

>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

```

The filtering condition is hard-coded. You could potentially parametrise it if you feel like later you might want to filter by a different number of unknown aa.
> 1 votes

# Answer
**Updated**: as pointed out, this will not work for a multifasta as asked for in your original question. However it would work fine for a single line fasta

---

You could accomplish this with a single line of bash code if you would like (this assumes your fasta is in gzip format, if not just use cat):

```
zcat prot.fasta.gz | grep -E -v -B1 '(^>.*|^(.*X.*){3}$)' | grep -E -v '^--$'  - | gzip > prot.filtered.fasta.gzip

```

The `-v` option of grep inverts the regex, and the `-B1` will print out the line occuring before the current context. The `-B1` introduces `--` between the matches, which is why the second grep is there. Then you can immediately pipe it back into a zipped file
> 0 votes

---
Tags: sequence-alignment, ngs, parsing
---


# Question
Title: Difference between PPI networks and Gene Co-expression networks
What exactly is the difference between protein-protein interaction (PPI) networks and gene co-expression networks? Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks, and so represent processed data. But could someone elaborate more on the differences?

# Answer
> 1 votes
> Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks

Not really, no. PPI stands for protein-protein interaction. In PPI networks, each node is a protein and each edge (line) connecting it to another node represents an interaction. A very common example is that an edge represents a physical interaction between these two proteins. This means that there have been experiments to indicate that these two proteins come into contact in the cell and interact.

For example, they might form part of a protein complex, so their interaction was seen using chromatin immunoprecipitation, where the two proteins were captured together. Or, they may have been seen together using FRET, a technique that lets you detect when two molecules are in very close proximity.

Alternatively, some interactions can be theoretical, where *in silico* modelling has suggested that these two proteins can interact, even though we don't have any direct evidence that they do so.

Finally, yes, it is also possible to make a PPI network based on gene co-expression data, under the assumption that if two genes are co-expressed, then their protein products will be present in the same tissue. However, that doesn't mean they interact in any way, so it isn't really a very common way of building PPI networks, as far as I know.

In co-expression networks, on the other hand, each node is a gene and not its protein product. And here, the edges indicate the coexpression. Often the level of correlation between the two genes' expression data.

The most important difference between these two kinds of network, apart from the obvious mentioned above, is that they are used in different ways. A gene coexpression network lets you identify genes that are regulated in similar ways. It is a tool you can use to investigate the regulatory network of the cell, to see what genes are "turned on" or "off" together.

The PPI network, on the other hand, is a way of understanding the result of this gene co-regulation. You are not looking at how the gene network is regulated any more, but how the products of the genes interact with one another. So it is a very different focus.

You can maybe think of this in terms of school children and their parents. You could make a network where each node represents a pair of parents, and then draw lines between them if their children were born in the same year. That would give you a general overview of when this group of people procreated. However, that two children were born on the same day doesn't necessarily imply that the children themselves will ever come into contact. This could be the equivalent of the gene co-expression network.

Then, make another network where each node is a child and draw a line between any children who know each other. This would be the equivalent of the PPI network. You would likely find some correlations between the two, with children who know each other being more likely to have been born at around the same time, but the two networks give you very different information.

# Answer
> 0 votes
Co-expression networks are based on gene expression data. Genes are compared using their expression profiles over multiple samples. A common distance measurement is the Pearson Correlation Coefficient, assuming that co-expressed genes (genes expressed simultaneously in similar conditions) will have a high correlation coefficient.

In its simple form, a co-expression network is done like this. An expression matrix containing n genes x p samples is used in input (from microarrays or RNA-seq) to calculate a distance matrix of size n x n. Then the distance matrix is thresholded to keep only the most significant relationship (eg correlation \>0.8). Gene pairs are connected if their correlation coefficient is above this threshold. Then all connections are visualized in a network, which can then be subjected to community detection algorithms to find more densely connected set of genes.

PPI are usually built on in vitro assays such as Yeast-2-Hybrid. But connections can also be inferred from literature data and expression data. You may have a look at the StringDB for example to see which kind of data can be used to connect proteins.

---
Tags: gene-expression, protein-protein-interaction
---


# Question
Title: Passing custom distance matrix to T_Coffee
I am passing a custom distance matrix to T\_Coffee (written in C), as per the docs, but can't seem to get it to work.

The documentation says to pass a matrix like this using BLAST format https://tcoffee.readthedocs.io/en/latest/tcoffee\_technical\_documentation.html#blast-format-recommended

```
# BLAST_MATRIX FORMAT
# ALPHABET=AGCT
A G C T
A 0 1 2 3
G 0 2 3 4
C 1 1 2 3
...

```

I've tried all combinations (with and without the leading # comments, and with and without the trailing) ..., but no matter what I keep getting this error:

> Parsing\_error when reading blast\_matrix doolittle: A C 0

Here, my distance matrix file is titled doolittle because it is the difference between the hydropathy indices of amino acids based on the Kyte-Doolittle scale * 10 to make them all integers. I've pasted the full matrix below.

I am calling the command from linux 18.08 ubuntu terminal (although I've tried most other combinations.)

```
t_coffee -in=rhomboids.fasta -matrix doolittle

```

Any help would be appreciated. Thank you!

doolittle

```
# BLAST_MATRIX FORMAT
# ALPHABET=ACDEFGHIKLMNPQRSTVWY
A C D E F G H I K L M N P Q R S T V W Y
A 0 7 53 53 10 22 50 27 57 20 1 53 34 53 58 26 25 24 27 31
C 7 0 60 60 3 29 57 20 64 13 6 60 41 60 65 33 32 17 34 38
D 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
E 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
F 10 3 63 63 0 32 60 17 67 10 9 63 44 63 68 36 35 14 37 41
G 22 29 31 31 32 0 28 49 35 42 23 31 12 31 36 4 3 46 5 9
H 50 57 3 3 60 28 0 77 7 70 51 3 16 3 8 24 25 74 23 19
I 27 20 80 80 17 49 77 0 84 7 26 80 61 80 85 53 52 3 54 58
K 57 64 4 4 67 35 7 84 0 77 58 4 23 4 1 31 32 81 30 26
L 20 13 73 73 10 42 70 7 77 0 19 73 54 73 78 46 45 4 47 51
M 1 6 54 54 9 23 51 26 58 19 0 54 35 54 59 27 26 23 28 32
N 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
P 34 41 19 19 44 12 16 61 23 54 35 19 0 19 24 8 9 58 7 3
Q 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
R 58 65 5 5 68 36 8 85 1 78 59 5 24 5 0 32 33 82 31 27
S 26 33 27 27 36 4 24 53 31 46 27 27 8 27 32 0 1 50 1 5
T 25 32 28 28 35 3 25 52 32 45 26 28 9 28 33 1 0 49 2 6
V 24 17 77 77 14 46 74 3 81 4 23 77 58 77 82 50 49 0 51 55
W 27 34 26 26 37 5 23 54 30 47 28 26 7 26 31 1 2 51 0 4
Y 31 38 22 22 41 9 19 58 26 51 32 22 3 22 27 5 6 55 4 0

```

# Answer
> 3 votes
So I didn't find an answer,but I did find a way around if it anyone else has this problem. It does involve modifying the source code, but even if you don't understand c (I barely do) its still easy to do.

1) Download source code from github https://github.com/cbcrg/tcoffee Source code is here, view the dependency list - you will need all these libraries to compile the code for the change to take effect.

2) Add your custom matrix to source code. Search `matrices.h`. In this version it is `t_coffee/lib/matrices.h` You will see how they are defining their default matrices. Just add your matrix with the same format and give it a name. For example:

```
int customMT[]={
  4,
  0,  5,
 -3, -2, 17,
  0,  5, -3,  9,
  0,  0,  1,  1,  6,
 -2, -3, -3, -5, -4, 10,
  0,  0, -4, -1, -2, -3,  8,
 -2, -2, -5, -2,  0, -3, -3, 14,
  0, -2, -2, -4, -3,  0, -1, -2,  6,
  0,  0, -3,  0,  2, -1, -1, -2, -2,  4,
 -1, -1,  0, -1, -1,  2, -2, -1,  2, -2,  4,
  1, -2, -2, -3, -1, -2, -2,  2,  1,  2,  2,  6,
  0,  4, -1,  1, -1, -1,  0, -1,  0,  0, -2,  0,  8,
 -1, -2, -3, -1,  1, -4, -1,  1, -3,  1, -3, -4, -3, 11,
  1, -1, -2, -1,  2, -3, -2,  0, -2,  0, -2, -1, -1,  0,  8,
 -1, -2, -2, -1, -1, -1, -2, -1, -3,  1, -2,  0, -2, -1,  3,  8,
  1,  0, -2,  0,  0, -1,  0, -1, -1,  0, -2, -2,  0, -1, -1, -1,  4,
  1,  0, -2, -1, -2, -2, -2, -2,  0, -1,  0,  0,  1,  0,  0, -3,  2,  5,
  1, -2, -2, -2, -3,  1, -3, -3,  4, -2,  1,  0, -2, -4, -3, -1, -1,  1,  5,
 -5, -5, -2, -4, -1,  1,  1, -5, -3, -2, -2, -3, -7, -3, -1,  0, -3, -5, -3, 20,
  0, -1, -2, -1, -1, -1, -1, -1,  0,  0,  0,  0,  0, -1,  0, -1,  0,  0,  0, -2, -1,
 -4, -3, -6, -1, -2,  3, -3,  0, -1, -1,  3, -1, -4, -2, -1,  0, -2, -1,  1,  5, -1,  9,
  0,  0,  0,  0,  5, -4, -2,  0, -3,  1, -1, -1, -1,  0,  4,  0, -1, -1, -3, -1,  0, -2,  4};

```

Then, search for file `tcoffee/lib/io_lib/io_func.c` Find the function

```
int ** read_matrices( char *mat_name_in)

```

You will see a long list of if else clauses, add the name of your custom matrix using the same pattern as they do.

```
/*Read Header Matrices*/
if (strm(mat_name, "pam250mt"))vector=pam250mt;
else if (strm(mat_name, "idmat"))vector=idmat;
else if (strm(mat_name, "dna_idmat"))vector=idmat;

/*add custom matrix here*/
else if (strm(mat_name, "customMT"))vector=customMT;

```

3) Save the files, move into tcoffee/src/. Rebuild the files using the make file as described on the github installation. It is here you need make, g++, perl, and whatever else on the list as dependencies

```
make t_coffee

```

4) Now, align your file! The only thing here is to use the t\_coffee of the source code you just edited, and not t\_coffee that might already be on your path. So change t\_coffee/src/t\_coffee to whatever the relative path to this folder is.

```
tcoffee/src/t_coffee rhomboids.fasta -matrix customMT

```

And that's it!

---
Tags: sequence-alignment, covid-19
---


# Question
Title: Get a single number representing the contact probability between a pair of genomic loci using HiC data
Hello to the experts in analyzing chromosome structure data in HiC format,

I have a very basic question. I have a HiC file (specifics are mentioned below), using which I would like to obtain contact probability between a pair of genomic loci. I obtained contact probabilities between the two loci at different resolutions available in the HiC file. However, rather than having a vector of contact probabilities, I would like to have a single number that appropriately represents the contact probability between the loci.

My guess is that there could be two ways to get a single number as the contact probability between two given genomic loci.

1. Should I take the highest contact probability? or,
2. Depending on the lengths of the genomics loci of interest, may be the resolution that is close to that could be selected as a single number that represents the contact probability between the loci. Eg. with this strategy, if genomic loci of interest are ~5kb in length the contact probability at 5kb resolution is selected as a single number. Then, the issue is that the length of the genomic loci of interest could be lower than the lowest length of resolution available in the HiC file.
or,
3. Contact probabilities at different resolutions could be aggregated. Then I wonder what would be the appropriate way to aggregate contact probabilities, e.g. by multiplication or median or something else.

Just to reiterate, above these are my guess and I am not sure whats the best way. So basically my question is: **how can I get a single number that would represent the contact probability between a pair of genomic loci using HiC data**?

---

Some specifics that might be useful in answering my question.

# Answer
> 2 votes
The fundamental unit of a Hi-C contact is the nearest restriction enzyme cut site. This is where each mate in a read pair from a Hi-C library comes from and is mapped to when converting sequencing data into a contact matrix.

Because the distance to the nearest cut site is variable at different genomic positions, and depth of coverage of any single locus is typically low, contacts are aggregated into bins and "normalized" to account for differences in cut site densities within those bins. This normalization transforms the raw contact matrix into something resembling a stochastic transition matrix, where you can think of the elements of that matrix as "what is the probability that at any moment in time bin $i$ is in contact with bin $j$".

So the probability you're looking for is the $(i, j)$th element of the normalized contact matrix, where bins $i$ and $j$ contain your regions of interest. If your regions straddle the bin boundaries, you may want to average those bins.

You get to pick the resolution of the contact matrix, but if you select a resolution that's too small, you run the risk of making your matrix too sparse and having no contact between those loci. If you pick one that's too large, your matrix elements aren't really informative for your question at hand.

Given your 500-5000 bp range and median of 1 kbp, it may make sense to pick a resolution between 1 kbp and 5 kbp and use that.

---
Tags: genome, hi-c, chromosome-capture, chromosomes
---


# Question
Title: Arlequin files not able to converge beyond 2000 steps for some .arp files
I am looking at 59 mitochondrial CO1 DNA sequences of length 563bp, but it is a haplotype file and it has 9 haplotypes within.

When I use .arp files in arlequin to run mismatch distribution it says that can’t converge beyond 2000 steps and also no graphs are produced. Can someone suggest what can I do to overcome this problem?

# Answer
> 0 votes
What you appear to have is a Bayesian MCMC for 2000 "steps" which isn't a lot and normally the term 'steps' is used for parsimony, so I'm a bit confused.

However, if a Bayesian analysis does not reach convergence, i.e. duplicate runs do not achieve the same results then the data space has not been solved so you get no result.

9 haplotypes does not sound much to me, it is possible there are too many sequences, it is also possible there is simply not enough variation (9 haplotypes ... how many alleles?).

Suggestions

* Arlequin should provide a default convergence value which you can lower to force a result
* Decrease the number of sequences but maintain the number of haplotypes
* Try their example data set

Do you get any associated file such as posterior probability, any output file with a minus sign in front of it?

---
Tags: statistics, sequence-analysis
---


# Question
Title: How to split FASTQ reads without re-running `fastq-dump`?
I downloaded a few FASTQs from the SRA that were paired-end 76 bp. When I look at the FASTQs, I get something like this:

```
@SRR7702334.1 1 length=152
AAGATNTTTTTCAGTATCAATCAGANGNNANANNCNNCTGCCCCAAATCCTGCAAAGAGCATGCAACGTCCTTAGAAGACAACGTCTTGCTCTGTCTCACCCAGGCTGGAGCATAGGGGCATGATCGATCATGACAAGGAGTTTACATTTTG

```

I forgot to add the `--split-3` option to my `fastq-dump` command. Is there an easy and fast way, without re-downloading the files, to split the reads in half and write them to two output FASTQs?

I tried the following using Biopython, and it'll work, but it's going to take a long time to finish.

from os.path import basename
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from gzip import open as gzopen
import argparse
from tqdm import tqdm

def split_read(read, pos=152 // 2):
    """
    Split read in half to produce 2 reads

    Parameters
    ----------
    read_id : Bio.SeqRecord.SeqRecord
        Read to be split
    """
    return {
        # R1/R2 as keys
        "R" + str(i + 1): SeqRecord(
            # first and second halves of sequence string
            seq=Seq(str(read.seq[(pos * i):(pos * (i + 1))]), IUPAC.ambiguous_dna),
            id=read.id,
            description="",
            letter_annotations={"phred_quality": read.letter_annotations["phred_quality"][(pos * i):(pos * (i + 1))]}
        ) for i in range(2)
    }

def split_fastq(fastq, prefix="split"):
    """
    Split reads from a FASTQ

    Parameters
    ----------
    fastq : str
        Path to input FASTQ file
    prefix : str
        Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)
    """
    # check for file format
    if fastq.endswith(".gz"):
        fq_handle = gzopen(fastq, "rt")
    else:
        fq_handle = open(fastq, "r")
    out_handles = {r: gzopen(prefix + "_" + r + ".fastq.gz", "wt") for r in ["R1", "R2"]}
    # load reads for random access
    records = SeqIO.parse(fq_handle, "fastq")
    for read in tqdm(records):
        # split read into each mate and fix header information
        mates = split_read(read)
        # write to each output file
        for m in mates:
            out_handles[m].write(mates[m].format("fastq"))
    fq_handle.close()
    for m in mates:
        out_handle[m].close()

if __name__ == "__main__":
    # parse command line arguments
    PARSER = argparse.ArgumentParser(description="Split a FASTQ into 2 FASTQs by splitting each read in half")
    PARSER.add_argument(
        "fq",
        type=str,
        help="Input FASTQ to split"
    )
    PARSER.add_argument(
        "-p",
        "--prefix",
        type=str,
        help="Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)",
        default="split"
    )
    ARGS = PARSER.parse_args()
    # run splitting function
    split_fastq(ARGS.fq, ARGS.prefix)

```
</code></pre>

# Answer
My comment above mentions the slow IO for the Python code (it was iterating at ~ 3600 reads/s, which, for a FASTQ with 500M reads would take ~1.6 d to complete), so I tried an alternative method that would cut down on the IO overhead. It's not super efficient, but I ended up using `zcat | awk` with a `switch-case` on the line number to stream the data.

```
# first pass to generate _R1 FASTQ
zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 1, 76); break; case 3: print "+"; break; case 0: print substr($1, 1, 76)} }' | gzip -nc > {output}_R1.fastq.gz

# second pass to generate _R2 FASTQ
zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 77, 76); break; case 3: print "+"; break; case 0: print substr($1, 77, 76)} }' | gzip -nc > {output}_R2.fastq.gz

```

I played around with trying to pipe each half of the combined read into its own file, but it was difficult to set that up and generalize it for the multiple FASTQs that I wanted to fix. I didn't feel like over-engineering something, and this was straightforward enough. It worked much quicker than I expected.

All in all this took ~4 h to complete, which means an iteration speed of ~ 70 000 reads/s, or about 20x faster than the Biopython solution. If I engineered a solution that would have only passed through the FASTQ once, it would've taken half the time, but in either case it was likely shorter than re-downloading the entire FASTQ.

But I will probably make use of `fasterq-dump` next time, and make sure I specify to split the read mates.

I love how how extensive the default biopython fastq parser library is, but the slow IO for large sequencing files makes it unusable in a lot of cases, unfortunately. I wonder if there's a good way to use C bindings to make this faster, or whether you just need to use a different package/language.
> 2 votes

---
Tags: fastq, sratoolkit, fastq-dump
---


# Question
Title: EGA remote download via snakemake fails due to EGA_CLIENT_ID specification
The example code provided by the Snakemake documentation for remote downloading from the EGA archive fails for me. I have defined all required environment variables, but an error comes up stating

```
{"error":"invalid_client","error_description":"Bad client credentials"}

```

I guess I am not defining the `EGA_CLIENT_ID` variable correctly? I am using Snakemake 5.16 and Python 3.7.2. Thanks a lot in advance for your advice.

Code:

```
import snakemake.remote.EGA as EGA

ega = EGA.RemoteProvider()
os.environ["EGA_USERNAME"] = "XXX"
os.environ["EGA_PASSWORD"] = "XXX"
os.environ["EGA_CLIENT_SECRET"] = "XXX"
os.environ["EGA_CLIENT_ID"] = "pyega3"
rule a:
  input:
      ega.remote("ega/EGAD00001002142/COLO_829_EPleasance_TGENPipe.bam.bai")
  output:
      "data/COLO_829BL_BCGSC_IlluminaPipe.bam.bai"
  shell:
      "cp {input} {output}"

```

# Answer
> 0 votes
Okay, I finally solved it. The reason was the `EGA_CLIENT_ID` variable. In the `pyega3.py` script, there is a definition for the `client_id` for pyega3 v3.0.44 ("f20cd2d3-682a-4568-a53e-4262ef54c8f4"). After setting the `EGA_CLIENT_ID` to this value login was successful.

```
os.environ["EGA_CLIENT_ID"] = "f20cd2d3-682a-4568-a53e-4262ef54c8f4"

```

I think this is a bug and I filed it on Snakemake's git repository - the official documentation just requires definition of the `EGA_USERNAME` and `EGA_PASSWORD` variables, and while definition of `EGA_CLIENT_SECRET` will probably need to be done by the user, the client ID should be pulled automatically from the pyega3 version installed locally.

Of note, after successful login, the example download from the documentation failed with the following:

```
Missing input files for rule a:
ega/EGAD00001002142/COLO_829_EPleasance_TGENPipe.bam.bai

```

However, downloading files that I have access to worked, indicating that the paths to the files from the documentation are wrong or are access-protected.

---
Tags: snakemake
---


# Question
Title: cellranger mkfastq with full path to --id flag
I have always used `cellranger mkfastq` to demultiplex 10x genomics runs manually, though recently the commands to do so have been incorporated into a script that should call the `cellranger` whenever it detects that a run is finished (and can be demultiplexed).

The main difference when calling manually compared to being called automatically by the script is **from where** the commands are being called.

So calling the commands manually **from the directory where the samplesheet is and where the outputs shall be saved** with:

```
/home/ngs/installed/cellranger-3.1.0/cellranger mkfastq --id=demux_10x --run=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX --samplesheet=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/20200507_NovaSeq_129_Chromium_Anne_Fred_Samplesheet.csv --qc

```

does the job.

but now cellranger is being called (automatically) from `/home/ngs`. The error message that I am getting from cellranger is:

```
/home/ngs/installed/cellranger-3.1.0/cellranger mkfastq --id=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/ --run=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX --samplesheet=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/20200507_NovaSeq_129_Chromium_Anne_Fred_Samplesheet.csv --qc
/home/ngs/installed/cellranger-3.1.0/cellranger-cs/3.1.0/bin
cellranger mkfastq (3.1.0)
Copyright (c) 2019 10x Genomics, Inc.  All rights reserved.
-------------------------------------------------------------------------------

/home/ngs/installed/cellranger-3.1.0/cellranger-cs/3.1.0/tenkit/bin/_mkfastq: line 531: __/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/.mro: No such file or directory
Martian Runtime - '3.1.0-v3.2.3'

Invalid name: /nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/ (only numbers, letters, dash, and underscore allowed)

```

According to the error message I cannot use slashes `/` in the `--id` flag.

Does anyone what I should pass to `--id` in order to get the outputs saved as `/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/`?

I already know that if I just pass `--id=demux_10x` the outputs will be saved at `/home/ngs/demux_10x`. (again, I am calling the script that calls `cellranger` from `/home/ngs/`)

Obviously I could change to the desired output directory i.e. `/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/` and do `--id=demux_10`, though I am trying to avoid that for now.

# Answer
> 1 votes
You should be able to use the `--output-dir` argument to specify a folder of your choice and use `--id=demux_10x` as before.

# Answer
> 0 votes
I suspect the error message is a bit exaggerated. Try removing the slash at the end of your --id

---
Tags: rna-seq, 10x-genomics, cellranger, demultiplexing
---


title: Getting weird plot with ggpubr package

text:
# Question
Title: Getting weird plot with ggpubr package
I have a data frame for two groups and 9 variables

```
> head(aa)
  treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
> class(aa)
[1] "data.frame"
> str(aa)
'data.frame':   1026 obs. of  3 variables:
 $ treatment: Factor w/ 2 levels "post","pre": 2 2 2 2 2 2 2 2 2 2 ...
 $ variable : Factor w/ 9 levels "Sig1","Sig2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ value    : num  0.117 0.189 0.219 0.128 0.171 ...
> tail(aa)
     treatment variable value
1021      post     Sig9     0
1022      post     Sig9     0
1023      post     Sig9     0
1024      post     Sig9     0
1025      post     Sig9     0
1026      post     Sig9     0
>

```

I am using this code

```
library(ggpubr)
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment")
> my_comparisons <- c("pre","post")
> p + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y = 50)
Warning message:
Computation failed in `stat_signif()`:
missing value where TRUE/FALSE needed
>

```

But I am getting this weird plot

While something like this is my goal

I really don't know where I am doing wrong

I used ggplot where I got something like this but I need statistics that is why I have stuck on this

I then used this code

```
> DF <- aa %>% pivot_longer(.,cols = c(Sig1,Sig2,Sig3, Sig4,Sig5,Sig6,Sig7,Sig8,Sig9), names_to = "var", values_to = "val")
> View(DF)
> library(ggplot2)
> library(ggpubr)
> ggplot(DF, aes(x = treatment, y = val, fill = treatment))+
+     geom_boxplot()+
+     stat_compare_means(comparisons = c("pre","post")+
+     stat_compare_means()+
+     facet_grid(var~., scales = "free")
+

```

But I got another weird plot

But changing the scale I got this but gives one statistics

While I need one stat for each variable between pre and post groups

This solution gives to statistics

```
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment",facet.by = "variable")
> p + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y =2)

```

When I edited my code

```
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment")
> p + stat_compare_means(aes(group = variable))
>

```

But again one stat

# Answer
The problem is the scale used: For the plot you called "weird" (first from the top), the scale is 50 and for the "ggplot only" (third from the top) the scale is 1.

You should play with the `stat_compare_means(label.y = 50)` bit, you can try setting the `label.y` parameter to 1.5 or 2.

EDIT BECAUSE OF THE ADDITIONAL REMARKS ON THE QUESTION:

Once again (and for sure the last time), I had to manually create example data as OP ignores repeated reminders on how to make a reproducible example with R.

```
library(data.table)
library(ggpubr)

text <- "
treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
1       post     Sig1 0.1173863
2       post     Sig1 0.1888243
3       post     Sig1 0.2191765
4       post     Sig1 0.1277735
5       post     Sig1 0.1711712
6       post     Sig1 0.0000000
1021      pre     Sig9     0
1022      pre     Sig9     0
1023      pre     Sig9     0
1024      pre     Sig9     0
1025      pre     Sig9     0
1026      pre     Sig9     0
1021      post     Sig9     1
1022      post     Sig9     1
1023      post     Sig9     1
1024      post     Sig9     1
1025      post     Sig9     1
1026      post     Sig9     1
"
my_data <-fread(text)

> my_data
      V1 treatment variable     value
 1:    1       pre     Sig1 0.1173863
 2:    2       pre     Sig1 0.1888243
 3:    3       pre     Sig1 0.2191765
 4:    4       pre     Sig1 0.1277735
 5:    5       pre     Sig1 0.1711712
 6:    6       pre     Sig1 0.0000000
 7:    1      post     Sig1 0.1173863
 8:    2      post     Sig1 0.1888243
 9:    3      post     Sig1 0.2191765
10:    4      post     Sig1 0.1277735
11:    5      post     Sig1 0.1711712
12:    6      post     Sig1 0.0000000
13: 1021       pre     Sig9 0.0000000
14: 1022       pre     Sig9 0.0000000
15: 1023       pre     Sig9 0.0000000
16: 1024       pre     Sig9 0.0000000
17: 1025       pre     Sig9 0.0000000
18: 1026       pre     Sig9 0.0000000
19: 1021      post     Sig9 1.0000000
20: 1022      post     Sig9 1.0000000
21: 1023      post     Sig9 1.0000000
22: 1024      post     Sig9 1.0000000
23: 1025      post     Sig9 1.0000000
24: 1026      post     Sig9 1.0000000
      V1 treatment variable     value

ggboxplot(my_data,
          x = "treatment",
          y = "value",
          color = "treatment",
          facet.by = "variable") +
  stat_compare_means(label = "p.format")

```

And this code generates:
> 3 votes

---
Tags: ggplot2
---


# Question
Title: Evaluating amateur bioinformatics results
I've might have created a new algorithm for finding patterns in DNA. The technique has not been used before, and it managed to find a 33-mer with corresponding complement that exists both in E. Coli and in Salmonella. I think it's pretty neat but I don't know how to evaluate the results. In addition, I've developed this algorithm after a bioinformatics course where only the DNA sequences of E. Coli and Salmonella were provided, so I would like to test it on other sequences and have it evaluated by people who know what they are doing.

How would one get in contact with the right people for evaluating amateur bioinformatics result?

# Answer
> 1 votes
You need to understand the biological significance of the 33-mer locus.

1. I'd blast your 33-mer on NCBI and see what the biological significane is of the locus. If might be a transposon/ promoter region.
2. Look at the translation of the protein and if its part of a protein the location within the protein structure.
3. Blast the entire locus including 5' and 3' regions of the 33-mer target
4. Carefully assess the results for chance, e.g. conserved protein and the synonymous mutations just happened to drop out (Poison distribution).

Once you gain a greater understanding, it is likely to have been identified before and what you might contribute is a different method to obtain a known result.

In summary, you need the next level part of the course about translating a genetic correlation into a biological result/ mechanism. This isn't research grade work at present in my opinion.

---
Tags: sequence-analysis, algorithms
---


# Question
Title: Effect of mutation in DNA sequence on transcription factor binding sites
How much does a single mutation/alteration of a nucleotide affect the presence of a transcription factor binding site (TFBS)? I am from computer science background(Obviously). I want to make a general assumption about the number of mutated bases in a TFBS to minimize the binding affinity?

# Answer
I haven't worked on TFBSs but I really, really doubt you will be able to do this. In mutations, like in real estate, it's location, location, location. It isn't a question of how many mutations but what mutations and where. TFBSs are usually described using position weight matrices, precisely because that's how you can show the importance of different changes in different positions.

Here are a couple of random Google image search results for "transcription factor binding sites logo":

(This was taken from here and you might want to read through it.)

(source: Front. Genet., 23 February 2016 | https://doi.org/10.3389/fgene.2016.00024)

The main take home messages from these images are:

1. Different positions have different importance and potential to disrupt binding. In the first image, I would expect that changing the T at position -3 to anything else would be pretty bad. However, changing an A to a T at -8 would likely be neutral.
2. Different TFs have very different binding sites and very different binding site lengths. Look at panel B of the second image. That one appears to be a very conservative site and likely any change would be deleterious. Conversely, changing any of the nucleotides of the 17bp spacer shown in panel A would likely not be important as long as you don't change the length.

Please note that I have not done more than glance at the sources I took these from, so my inferences may well be wrong, I only use them to illustrate the general point: you cannot really put a number to this. You will always need to look at each specific change and each specific position. In some cases, a single change can disrupt binding while in others multiple changes might have no effect at all.

---

You might want to look into TFBS databases like JASPAR or TRANSFAC.
> 2 votes

# Answer
Every motif has a score. The score is calculated from the number of times each base was observed in that motif in the genome (I'm not sure how exactly these are calculated). Higher numbers indicate that this motif was found more frequently, which we take as a proxy for it providing better binding.

The Ensembl VEP can give you the score change for a mutated motif in the output, so this may be the easiest way to analyse this.
> 1 votes

# Answer
The motifbreakR package for R I believe does something like what you want to do;

https://www.rdocumentation.org/packages/motifbreakR/versions/1.2.2/topics/motifbreakR

If you're looking for a general rule as to what kind of mutations/where in the motif mutations will significantly affect motif binding, I suggest generating a lot of sequences with simulated SNP's (single nucleotide polymorphisms) and seeing how these score against the particular motif. Then, you will have a large(ish) database of sequences with one base pair difference from the motif and their corresponding score. This will allow you to judge what kinds of mutations affect the score the most.

Furthermore, you could somehow make a distribution of scores versus mutations and generate an empirical p-value to determine statistical significance (for instance, such and such mutation caused a score that was in the lowest 5% of scores for all mutations). You would know better than me how to do this, since your background is in CS and my PhD subject is genetics.

To make an assumption about the number of bases that need to be mutated, do this for one, two, three, etc. simulated SNP's.
> 0 votes

---
Tags: fasta, sequence-analysis, phylogenetics
---


# Question
Title: Extracting specific sections from FASTQ determined by positions
Given a FASTQ file, I'd like to generate a new FASTQ file including only certain subsections of the original sequences specified by their positions.

For example, say I want to extract the nucleotides at position **1, 3, 5-7** from each sequence:

Input:

```
@id1
AACCGGTCC
+
123456789
@id2
TACCGGCCC
+
123456789

```

Output:

```
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

```

Note, the phred-scores are constructed such that they correspond to the position for easier reading. Note2, the phred-score subsections should be included accordingly (i.e. the same positions)

Ideally, a command should work this way:

```
subselect 1,3,5-7 input.fq > output.fq

```

Having looked through several poplar software, I couldn't find anything working out of the box. I was thinking to use either awk or even cutadapt (systematically trimming from the 5\` and pasting) but maybe out of an box solution does already exist which evaded me.

# Answer
> 1 votes
For the specific regions you mentioned, you could do this in `awk`. It won't be fast for a large file, but it will work:

```
$ awk 'NR%4==2 || NR%4==0{$0=substr($0,1,1)substr($0,3,1)substr($0,5,3)}1 ' file.fq
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

```

`$0` is the current line. The script above will replace the line with three substrings of itself, one starting at position `1` with a length of one, one starting at position 3 with a length of one and the last starting at plsition 5 with a length of 3. The `NR%4==2 || NR%4==0` ensures that this substitution only occurs on lines whose line number modulo 4 is either 2 or 0, so every 2nd and 4th line of the file.

**IMPORTANT**: this assumes you *always* have only one line of sequence in each section. The fastq format allows multiple lines of sequence (and quality scores). So your data might not conform to this. However, for regular short-read data, it should be fine.

---

And here's a slightly more sophisticated version in Perl, with some rudimentary error checking and which can take regions in the format you asked for (`1,3,5-7`):

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

if ($#ARGV != 1) {
  die "Need exactly 2 arguments.\n";
}
my $file = $ARGV[0];

my @regions = split(/,/,$ARGV[1]);
my (@startPositions, @lengths);
for my $arg (0..$#regions) {
  if ($regions[$arg] =~ /^\d+$/) {
    push @startPositions, $regions[$arg];
  }
  elsif ($regions[$arg] =~ /^(\d+)-(\d+)$/) {
    for my $i ($1..$2) {
        push @startPositions, $i;
    }
  }
}

open(my $inputFileHandle, '<', $file) or
    die "Failed to open $file for reading: $!\n";

while (<$inputFileHandle>) {
  if ($. % 4 == 2 || $. % 4 == 0) {
    my $line = "";
    foreach my $start (@startPositions) {
      $line .= substr($_, $start-1, 1)
    }
    print "$line\n"
  }
  else {
    print
  }
  ## If we find a line that consists only of '+' but its line number
  ## modulo 4 isn't 3, that means that there was at least one multi-line
  ## sequence that we cannot handle correctly.
  if (/^\+$/ && $. % 4 != 3) {
    die "Line $. is just a '+', multi-line sequence found in file.\n"
  }
}

```

Save that as `foo.pl` and then run like this:

```
perl foo.pl file.fq 1,3,5-7

```

# Answer
> -1 votes
This looks like a job for samtools fqidx (http://www.htslib.org/doc/samtools-fqidx.html), assuming that your file is not enormous... read the warning in the link above.

usage would be

`samtools fqidx yourfile.fq <region>`

where region is specified using the fastq record identifier and coordinates.

---
Tags: fastq, awk
---


# Question
Title: Differential gene expression bias due to effect of an individual sample
I am analysing a human single cell RNA seq experiment, where we have 4 groups, four samples each. Data has been analysed using Seurat, with the canonical workflow. I have tried DE using various methods (Wilcox, MAST, LR), and for each of them the top genes (biggest fold change, lowest adjusted p) are genes over expressed by specific samples, e.g.:

Looking at the average expression for these genes per each subject, and plotting on a Boxplot, it does not seem to be a genuine difference (e.g. for the gene above):

Besides the subject, there is no other parameter (e.g. age, sex) that is biased across groups and that shows up as differentially represented in the PCA plot.

Is there a DE method that takes into account this behaviour and corrects for it? How can these results be presented? Besides increasing number of subjects (not really feasible), is there a solution to analyse these data?

Thanks a lot!

# Answer
Leave-one-out cross validation is an iterative method that leaves out one sample until each sample has been left out once. You can then identify the intersection of the differentially expressed genes from the sample space of repeated analyses to produce a set of genes that are consistently differentially expressed independent of individual sample biases. more info here: http://efavdb.com/leave-one-out-cross-validation/
> 3 votes

# Answer
One potential problem here is that, as far as I can see, you're not taking the hierarchical structure of your scRNA-seq data into account. I.e. the fact that cells from the same patient / sample will be more similar than cells from different samples, even within the same group. This is also the problem addressed in the paper mentioned by @alan. Also, it usually makes more sense to do this kind of group vs. group comparison **for each subpopulation (cluster/cell type) separately** (assuming you have multiple subpopulations in your data and did some kind of clustering).

One way of doing this is by using a **pseudobulk-level approach**, where you basically aggregate (by either averaging or summing) your data over each sample, for each cluster separately. So you end up with one counts matrix per subpopulation with one row for each gene and one column for each **sample**. The sample columns are then the sums (or averages) of all cells belonging to that cluster-sample combination.

You can then use bulk-RNAseq methods such as **edgeR** or **DESeq2** on the aggregated data to test differences between groups of samples for each cluster. The muscat package implements all these steps and has quite extensive documentation. So you could check that and test it out on your data.

The paper describing the method of *muscat* can be found here: https://www.biorxiv.org/content/10.1101/713412v2
> 2 votes

---
Tags: rna-seq, scrnaseq, seurat, differential-expression
---


# Question
Title: Should biological replicates be the most similar pairs in an RNAseq experiment?
The attached figure (from deseq2) shows the sample to sample distances of a RNAseq experiment with 4 conditions (A,B,C,D) at 2 time points (0h,4h) with 2 biological replicates.

I am a little bit puzzled that the biological replicates are not always the most similar pairs. In a perfect world biological replicates should be identical but I assumed that at least they should be very similar.

Is this something I have to worry about, or is this more or less normal?

# Answer
In summary you basically have noisy data and to simply answer your question it depends on the genetic diversity of the samples you are dealing with. Note, from my recollection this was confusing to start with because the *"replicates"* column was labeled as *"treatment"*, but I could be wrong.

---

Revised answer, firstly there is zero stationarity - thus its a straight time-series.

**Summary** As a time series

* The statistical noise between A and D is lower than the difference in phenotypic expression between each and all othe replicates in the analysis ABCD - FOR DAY 0 VERSUS DAY 4, so they cluster and are distinct.
* The statistical noise between B and C is greater than the phenotypic variation between their expression profiles for DAY 0 to DAY 4. I suspect they are genetically very similar. Thats the conclusion, you can spice it up with different analysis, but I'd leave it because might find a lack of statistical significance in your core observations (A and D).
* Within the time-points the data is very noisy and doesn't appear to show a distinct signal. Thats okay, this is essentially a time series.
* Sample D is distinct in all scenarios (which is good)

**Result** To reiterate and expand a bit, the result you have is that sample A and D biological replicates are consisent within the analysis you are performing, i.e. day 0 verses day 4 (below). Thus they cluster together and this is proof of distinct phenotypes above the statistical noise (A with A and D with D). An intra-time point dendrogram would have a different result. The replicates samples B and C are not consistent and the interpretation is that the phenotypic expression profiles are so similar they are essentially identical within the boundaries of statistical noise.

**Analysis** The dendrogram is combines day 0 to day 4 analysis with day 0 to day 0 analysis and day 4 to day 4 analysis. HOWEVER, the clustering algorithm will correct for this and separate based on the largest distance (day 0 versus day 4). The signal that you get in your dendrogram is the largest signal in the data set which is day 0 vs day 4, and this is how UPGMA clustering works (I probably have previous posts on the statistical algorithm of UPGMA). So actually that is fine.

**Advice** With data this noisy I'd go no further in attempting to deploy critical probablity (a *p* value) because this data could be the result of a complex distribution. This would require a statistician. Personally I'd upload the PCA as a separate post, but its just my opinion, because UPGMA versus eigenvalue based PCA are statistically quite different approaches.
> 1 votes

---
Tags: rna-seq, deseq2
---


# Question
Title: How to understand and analyse RNA-seq data (for a beginner)?
I am trying to understand expression of a certain protein across Pseudomonas species. I downloaded an SRA file from NCBI and converted it to a fastaq file. I am not able to understand how to interpret and analyse this file. I want to get the mRNA sequence and see which genes are being expressed and find the mRNA of interest.

How do I proceed with getting the mRNA sequences and then for mapping them to genome?

# Answer
SRA files are raw sequencing data that have been converted to the NCBI storage format. These can be converted to fastq files, e.g. with fastq-dump from the SRAtoolkit.

Beyond that I suggest you start with the Hitchhiker's GUide to RNA-seq and the DESeq2 workflow to get a background. Be sure to extensively browse the web for vignettes, guides and tutorials. YOu will find helpful resources to get started.
> 1 votes

---
Tags: rna-seq, rna, transcriptome, gene-expression
---


# Question
Title: Tools for creating gene coexpression networks
What are some good tools (in Python if possible) for creating gene co-expression networks from RNA seq data?

# Answer
After looking around, it seems like the Camoco library is suitable.
> 1 votes

# Answer
it depends on the size of your data. If you have enough RAM available, then you may directly compute a correlation-based distance matrix using the *cor()* function (with method="spearman" or method="pearson"). But you may also give a try to Mutual Information distance using the **Parmigene** (https://cran.r-project.org/web/packages/parmigene/index.html) R package or partial correlations with the **Corpcor** (https://cran.r-project.org/web/packages/corpcor/index.html) R package. If you have a large expression matrix, remember to use the excellent **data.table** package and its *fread()* function.
> 1 votes

---
Tags: python, gene-expression
---


# Question
Title: local BLAST error: BLAST Database error: Error: Not a valid version 4 database
I want to use a local `blastn` command to BLAST a multifasta file ~700 sequences using following command: `blastn -db nt -query fasta_all.fasta -num_alignments 2 -out fasta_blasted.txt` but I receive error :`BLAST Database error: Error: Not a valid version 4 database.`

I use a local `nt` database which I downloaded today (11th of May 2020) by running `update_blastdb.pl --decompress nt [*]`. I downloaded this but I did not do any post-processing and just the "raw" download files are in the folder (see image below which shows only a part of the files)

The `blast+` I downloaded and installed is `ncbi-blast-2.10.0+-x64-linux.tar.gz` from https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

So, both database and `blastn` version should be the latest. It probably has something to do with either the version of the database or the software, but as it both are the latest version, I find it strange that these would not be compatible. Any ideas what causes this error?

# Answer
> 2 votes
first I recommend you to check your downloaded data with md5 to be sure there was no problem with that step. How many files did you download?

Also I was curious if you had another BLAST+ version installed on your computer. Be sure that blastn is called is a version corresponding to the database you are using.

---
Tags: database, public-databases, blastn
---


# Question
Title: Is it important to filter out poor quality cells before performing an integration analysis on single cell RNA sequencing data?
In order to perform an integration analysis of single cell RNA seq data, is it important to check the percentage expression of mitochondrial genes of cells as well as the feature counts to exclude dying or dead cells as well as multiplets from further analyses?

I ask because the standard Seurat work flow to perform integration as outlined in Integration and Label Transfer does not mention filtering steps prior to integration.

In this publication, Eleanor Alenzi et al. 2018 also did a canonical correlation analysis with Seurat but did not mention steps to filter out poor quality cells before doing the integrated analysis. This second publication does same. Both publications do not mention if their visualisations or differential gene expression analyses are done on integrated (discouraged) or raw uncorrected RNA data (recommended).

Following the styles of the publications linked above (many more do same), is it also efficient to perform integration without filtering out poor quality cells?

# Answer
> 4 votes
It is absolutely necessary to remove low quality cells: In the case of `CCA` (and this applies to other "integration" or "data alignment" methods as well), one would need to use "anchors", basically same type/state/kind of cells from the samples to be "integrated" and are used to "align" the different samples. In a scenario where you have a set of cells whose distinguished feature is high mitochondrial content due to apoptosis irrespective of their origin (epithelial, T cell, ...) and such cells are picked up as so called anchors, the whole integration effort will be off. Even when not using integration methods, bad quality cells quite often form their own clusters (highly skewed biology) not contributing anything to analyses and are often removed at the downstream steps anyway.

And when it comes to published work, unfortunately it comes down to the backgrounds of the reviewers. I tend to think that there are way too many single cell studies than experts with required bioinformatics and statistics so that a thorough review can be conducted. And lack of a paragraph on QC does not mean no QC has been performed anyway.

---
Tags: scrnaseq, seurat
---


# Question
Title: Identify differentially covered genes only between two samples
I have a question about finding differentially covered regions (coverage represents methylation level which goes from 0 to several thousands). I'm using enrichment based method which can be summarized with coverage per gene:

```
data <- matrix(sample(80), 20)
# Genes in rows
rownames(data) <- letters[1:20]
colnames(data) <- c("group_A_tr1", "group_A_tr2", "group_B_tr1", "group_B_tr2")

```

In `data` matrix each row represents a gene and each column represents a sample. There are two sample groups (*A* and *B*) with two technical replicates per each group. **Problem is that we do not have any biological replicates.**

My goal is to identify genes that are differentially methylated between two groups. I know that `limma`, `edgeR`, `DESeq2` can be used in analysis like this, however I don't have enough samples. Basically I'll need to compare only two columns (after averaging technical replicates).

What method would be appropriate to work with data like this? Is it possible to treat technical replicates as biological ones?

# Answer
Not really an answer but an extended comment... and most likely something you don't like to hear

I guess by technical replicates, it means taking the same biological sample and making 2 methylation libraries. If this is used as replicates in deseq2 or edger, the variance you are estimating is the technical variation that comes with preparing the library, bisulphite conversion etc..

Often the biological variance will be much larger. That is if you take 2 independent samples of group A, their variation in methylation is much larger than the technical variation you see from handling the samples. If the technical replicates are treated as biological replicates and differential analysis is performed, it will be seriously misleading.

I think both edgeR and DESeq2 doesn't estimate dispersion with 1 v 1 (which is right).

I would suggest maybe the following:

1. Run a fisher.test at every gene (methylated / unmethylated), get the odds ratio and p-value, rank this, but explain that the p-value is only for ranking, and will underestimate the variation.
2. A bit more refined, look for other studies that might be close to yours, get an estimate of the dispersion in their study and use a negative binomial distribution to again rank genes that show the biggest difference
> 4 votes

# Answer
You may run edgeR for methylation analysis without replicates (https://f1000research.com/articles/6-2055). But I also recommend you to have a look at the R DSS package (http://bioconductor.org/packages/release/bioc/html/DSS.html). It has a smoothing step which allows to modelize the biological variability by taking information from neighboring sites.
> 1 votes

---
Tags: bioconductor, deseq2, limma, edger, methylation
---


# Question
Title: Should you filter GWAS hits with high standard error?
I'm trying to figure out if I should be filtering out GWAS hits that have high standard error and I'm not quite sure what to do. It seems like it might not matter, because the standard error is used to calculate the t-statistic, which is then used to calculate the p-value. So in a way it's already built in. But reporting SNPs that have very high standard error doesn't seem quite right. What's the best way to handle this?

# Answer
I am curious about the relationship between your p-value (or effect sizes) and standard error. I would expect the significant signals to have smaller stand error compared to the non-significant, background signals. If this is the case, there is no need to do filter.
> 1 votes

# Answer
I think the effect size is more of an issue than the standard error.

If the standard error suggests that the effect direction could change sign, then it might be a good idea to filter it out. Otherwise, if the effect size remains large (or at least positive) even after accounting for a large error, then \[to me\] that would attach more weight to the likelihood that the association is valid.

I did a very quick analysis of the results of one UK Biobank study on Twitter, where I used the statistic of "the least extreme value in a 95% confidence interval" in a Manhattan plot instead of the more commonly-used p-value:

https://twitter.com/gringene\_bio/status/1207617723586371584
> 0 votes

---
Tags: statistics, gwas, genetics
---


# Question
Title: Importing a WGCNA co-expression network into a networkX graph in Python
Is there a way to generate a WGCNA co-expression graph in R, and then to import that graph into Python such that I can represent it using a networkX object?

# Answer
> 2 votes
You can use `exportNetworkToCytoscape` or `exportNetworkToVisANT` to generate and egde list file with the two nodes at the ends and weight for each edge. Hopefully there's a way to import that into networkX.

---
Tags: r, python, networks, wgcna
---


# Question
Title: When should I apply a multiple comparison adjustment?
I am performing a RQL analysis ("http://www.esapubs.org/archive/ecol/E095/002/suppl-1.pdf"). This analysis meassures the relationship between different traits of three different matrix showing if there is a positive or negative relationship. This RQL analysis make a fourth-corner method.

Well, this method have the option to make a p-adjust method to solve multiple comparison problem. However, my variables are independent and I belive that I should not apply to my dataset. Am I right? Or always that two or more matrix are compared I should apply a multiple comparison adjustment?

Thanks in advance.

# Answer
The p-value adjustments corrects the inflated number of significant results when performing many significance tests.
If you have a test with 5 % false positive rate and but you run it several times the chances of getting a false positive results in any of these tests is far greater than 5 %. For example repeating the test with 5 % false positive rate 20 times means that there is a 60 % chance of at least one false positive result.

You can see this in the tutorial you have linked to in your question, where the number of significant differences goes down from 26 without correction to 18 with correction.

Therefore, if you are performing multiple comparisons you should do a correction. The problem of inflated p-values will occur when you are performing many tests even if the tested variables are independent.
> 3 votes

---
Tags: statistics
---


# Question
Title: VCF-merge fails due to tabix not producing .tbi files
I'm trying to use vcf-merge to combine 2 exome capture vcf files (~250K and ~330K in size) before trying it on all 96 samples. I'd appreciate any advice on the best way to do that! I've detailed what I've tried below. My issue seems to be with using tabix to convert the files to .tbi format.

Step 1: BGZIP

So far, I've zipped the files without issue:

```
bgzip sample1.vcf
bgzip sample2.vcf

```

which produces:

```
sample1.vcf.gz
sample2.vcf.gz

```

Step 2: TABIX

When I try to use this command:

```
tabix -h -p vcf sample1.vcf.gz

```

the stderr is:

> Region 536999277..536999278 cannot be stored in a tbi index. Try using a csi index with min\_shift = 14, n\_lvls \>= 6. tbx\_index\_build failed: sample1.vcf.gz

Using the -C option which works:

```
tabix -C -h -p vcf sample1.vcf.gz
tabix -C -h -p vcf sample2.vcf.gz

```

which produces: sample1.vcf.gz.csi sample2.vcf.gz.csi

Step 3: VCF-MERGE

When I use this command:

```
vcf-merge sample1.vcf.gz.csi sample2.vcf.gz.csi > out.vcf.gz

```

the merge fails, and I get this stderr:

> Broken VCF header, no column names? at /usr/share/perl5/Vcf.pm line 172, \<**ANONIO**\> line 1. Vcf::throw(Vcf4\_2=HASH(0x5645a760bc38), "Broken VCF header, no column names?") called at /usr/share/perl5/Vcf.pm line 867 VcfReader::\_read\_column\_names(Vcf4\_2=HASH(0x5645a760bc38)) called at /usr/share/perl5/Vcf.pm line 602 VcfReader::parse\_header(Vcf4\_2=HASH(0x5645a760bc38)) called at /usr/bin/vcf-merge line 183 main::init\_cols(HASH(0x5645a761f438), Vcf4\_2=HASH(0x5645a760b248)) called at /usr/bin/vcf-merge line 279 main::merge\_vcf\_files(HASH(0x5645a761f438)) called at /usr/bin/vcf-merge line 12

If I exclude the tabix files, the merge still fails and the stderr says:

> The column names not tab-separated? Could not load .tbi index of sample1.vcf.gz. The command exited with an error. Is the file tabix indexed?

Any suggestions on how I can correctly merge these files?

# Answer
You should go with the CSI index as it seems like the BAI indexes cannot handle chromosomes longer than 2^29.

Your problem is that you're trying to merge the index files using `vcf-merge`. You should merge the VCF files, not the `.csi` index files.
> 2 votes

---
Tags: vcf, linux, vcftools, tabix
---


# Question
Title: how to change the UMAP use in the dimplot and feature plot
I would like to know how to change the UMAP used in `Dimplot` and `FeaturePlot` from Seurat: how we can get the x-axis and the y-axis like UMAP-1 and UMAP-2 if I want to use UMAP-4 and UMAP-5. How to modify the code? It is easy to change the PC by using `DimPlot(object = pbmc_small, dims = c(4, 5), reduction = "PCA")` but if I changed to `reduction = "UMAP"`, I got the error

```
"Error in Embeddings(object = object[[reduction]])[cells, dims] : subscript out of bounds Calls: DimPlot Execution halted".

```

Does anyone have any idea how to solve this problem?

# Answer
Try this:

```
library(Seurat)
pbmc_small <- RunUMAP(pbmc_small, dims = 1:10, reduction = "pca", n.components = 5)

names(pbmc_small[["umap"]])
#> [1] "UMAP_1" "UMAP_2" "UMAP_3" "UMAP_4" "UMAP_5"

# plot
DimPlot(pbmc_small, reduction = "umap", dims = c(4, 5))

```

This should work with Seurat version 3.1.5

```
packageVersion("Seurat")
#> [1] '3.1.5'

```
> 1 votes

# Answer
Change `n.components` / `max.dims` (depending on your seurat version) to 5 when you call `runUMAP` and you should have dim 4 and 5 available.
> 0 votes

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


title: RepeatMasker annotation

text:
# Question
Title: RepeatMasker annotation
I collected my own denovo library using many different tools Pipeline (structure-based ) and I updated all the headers corresponding to RepeatMasker format, to obtain Repeat coverage in the genome and .tbl annotation.

when I run repeatmasker the annotation results assume to be Homosapien, where Homo annotation table ( .tbl) is limited as shown below Table (1)

Actually my denovo repeat library consists of more than what is in the table above, example, rolling circles, Copia, Gypsy etc..

```
RepeatMasker -Lib myOwnLibrary.fa Mygenome.fa

```

Table (1)

My question is how can I run RepeatMasker with my own library using -lib Flag and result like the table (2) below,

the table below assume to be vertebrates genome, this is my aim.

```
RepeatMasker -Species Vertebrates Mygenome.fa

```

Table (2)

Regards

# Answer
> 1 votes
After long years of just conceptual knowledge of what the software does, I was dealing with `RepearMasker` this week. So, I am not entirely sure I am correct, but this is how I understand it.

The two flags you explored serve for different scenarios:

* `-species` \- if the taxa/species is present in the repeat database (`repBase` or `Dfam`)
* `-lib` \- if you want to give repeatModeler something different.

In your case it seems like your library is completely missing all LINEs, SINEs and a few others. What you can do is to get from the database all repetitions of the corresponding taxa and catinate it with the repeats you are interested in (the one you supplied in the `-lib` argument in the first command).

```
util/queryRepeatDatabase.pl -species Vertebrates > flavoured_vertebrate_repetitions.lib
cat myOwnLibrary.fa >> flavoured_vertebrate_repetitions.lib
RepeatMasker -lib flavoured_vertebrate_repetitions.lib Mygenome.fa

```

This should produce annotation similar to the bottom case, but also using the repetitions you manually added.

P.S. `util/queryRepeatDatabase.pl` is a script to query the database of repeats and it is shipped together with `repeatModeler`.

-- edit --

The TE classification in the summary file of repeatMasker is based on headers of the repetitive elements in the fasta file you specify in to `-lib` argument. So if you want to get the second-table-like, you need to somehow reclassify the repeat library to contan repBase-like classification of TEs. The explanation of lib argument from the manual:

> The recommended format for IDs in a custom library is:
>
> \>repeatname#class/subclass or simply
>
> \>repeatname#class
>
> In this format, the data will be processed (overlapping repeats are merged etc), alternative output (.ace or .gff) can be created and an overview .tbl file will be created. Classes that will be displayed in the .tbl file are 'SINE', 'LINE', 'LTR', 'DNA', 'Satellite', anything with 'RNA' in it, 'Simple\_repeat', and 'Other' or 'Unknown' (the latter defaults when class is missing). Subclasses are plentiful. They are not all tabulated in the .tbl file or necessarily spelled identically as in the repeat files, so check the RepeatMasker.embl file for names that can be parsed into the .tbl file.

# Answer
> 2 votes
RepeatMasker comes with a buildSummary.pl script.

You can input the .out file from your RepeatMasker run and there is a flag to enter your species

e.g.

```
buildSummary.pl -species Vertebrates rmrun.out > mysummary.tbl

```

In the mysummary.tbl, the summary will not assume your run to be *Homo* (default).

---
Tags: sequence-annotation, repeat-elements, repeatmasker, repeat
---


# Question
Title: umap and Louvain clustering on normalized data
I know that scaled data must be used for PCA for example as it is based on variance maximization.

However I'm wondering if it's the same case for UMAP ? If the data are single-cell RNA seq, after normalization can we do UMAP and Louvain clustering ?

I tried both and get similar results, however the Louvain clustering seems to be more adequate on normalized data than on scaled data.

By adequate I mean the clusters are the same but some are split into two, which makes sens looking at other results (like transcription factor analysis).

# Answer
> 3 votes
> If the data are single-cell RNA seq, after normalization can we do UMAP and Louvain clustering ?

Typically people run PCA, UMAP and Louvain clustering on the normalised and log-transformed expression counts (but do marker gene and differential expression analysis on the non-normalised values). However, this remains controversial. I recommend reading "Current best practices in single‐cell RNA‐seq analysis: a tutorial"\- it's a bit outdated, but still a good reference.

> By adequate I mean the clusters are the same but some are split into two, which makes sens looking at other results (like transcription factor analysis).

The number of (or size of) clusters should not be used to interpret how good the clustering methodology is, since it is generally dependant on the input parameters, e.g. a higher resolution parameter leads to more clusters. Instead, display the clusters on the UMAP and see if they are continuous 'visually' and check what genes are differentially expressed between clusters and see if they make sense 'biologically'.

# Answer
> 2 votes
There is discussion if scaling the data (making the data range for all genes the same) is something you need for single cell data. The argument is that the expression differences between genes themselves are informative. Normalization on the other hand is always necessary.
The big single cell pipelines like `Seurat` or `Monocle` use both normalization and scaling as standard.

If you would like some of the clusters you get split into subclusters you could either increase the resolution parameter that the clustering functions usually have. Or you can take the cells from a cluster, then subcluster them and add the subclusters to the original clustering results.

For yor next question it would be great if you could provide the plots of the data because this way it is easier to understand what the problem is.

# Answer
> 0 votes
If you working on known cell culture, I recommend using classification approach instead of clustering. One library that I used to publish my paper is SingleR. I checked the expression x cell type matrix by both approaches (Seurat Clustering and SingleR), the classification approach make much more sense.

---
Tags: scrnaseq, single-cell, umap
---


# Question
Title: Standalone Blast+; Automating searches, Blastn formatting, interpretation
I am struggling with using NCBI's standalone blast, particularly the blastn feature.

I have a large number of nucleotide sequences that I want to use blastn to ID. The goal is to retrieve information from the blast program and eventually import this information to excel. I want to know the capabilities of, or the easiest way, to format the output. I am running a remote search to NCBI's nucleotide database and want to filter and organize the information retrieved as much as possible. The ultimate goal is automate a way to run a remote search to NCBI's nucleotide database, and to work with the results in excel.

I am extremely new to this and struggling to find the best path. From what I can gather, many prefer using external programs, scripts, parsers, interpreters etc to organize their results. I am okay with doing this as well, but would still like to format the output of my blastn search as much as possible beforehand. It is returning a lot of data to me that I don't need.

Can I automate this in excel, and just use VBA/macros to further organize my data? It seems some people prefer to just import the output text file to a spreadsheet and work from there. I have also seen options like BioEdit, eutilities, APIs... I'm not sure which route I should invest my time into figuring out.

Thanks so much for lending your time and expertise!

# Answer
You just have to change the output format using the `-outfmt` parameter. To have an output in a tabular format just use `-outfmt 6`. The resulting table have 12 columns containing the following informations *queryacc subjectacc pident length mismatch gapopen querystart queryend subjectstart subjectend evalue bitscore*. You may restrict the number of hits per query using `-max_target_seqs 5` for example.
> 2 votes

---
Tags: blast, blastn
---


# Question
Title: Best practice to perform a functional analysis on a set of VCFs from a cohort of patients
I am confused about how to perform a functional analysis after a variant calling on a cohort of patients. I have anotated all vcfs and filter them in order to get only those deleterious variants.

I was thinking on merging those vcf: I tried to merge them with `bcftools merge` and also using `VariantAnnotation` and transforming each sample vcf into a dataframe containing the variant lines to then row bind all those dataframes together.

As an example, Ill show some fields of teh first row of a dataframe:

```
CHROM   START REF ALT FILTER DP  GENEID  SIFT
  17     1000   T   A   PASS   62   TP53    0

```

Then, I would extract the genes affected by those variants and perform a gene set enrichment analysis, a pathway analysis similar to those performed after an RNA-seq analysis (with `enrichR`, for example). But I am confused about the viability of this steps. For example, I would like that, during the pathways analysis, the program I use not just reads the list of mutated genes in my cohort, but also takes into account that, if a gene is mutated in more than one sample, it should wheight more in the analysis than a gene which appears in just one sample.

# Answer
As a simple approach you may want to try the "fuzzy" list implemented in EnrichR itself. There is some information about this in the Help Center under the "What is a gene set" heading.

This way you can code the number of data sets a gene appears as a value between 0 and 1. In the Help Center this is used to convert fold changes or p-values to a membership score but I don't see why this would not also work for the number of data sets a SNP is observed in.

Since it is relatively easy to implement this can give you an first impression of what to expect. It then may be combined with analyzing each data set independently and then comparing the overlap.
> 1 votes

---
Tags: variant-calling
---


# Question
Title: Hypergeometic test for the overlap of two set of genomic intervals (e.g. from Chip-seq data)
Suppose I have done two Chip-seq experiments. Now after the bias correction step such as MACS2, I get two set of genomic intervals (or peaks, as what is usually called) that corresponds to each of the experiment. I want to test the significance of the overlap of these two set genomic intervals that corresponds to each experiment.

The common test employed seemed to me is the Hypergeometic test. There's even a package specifically aim for this task: ChiPpeakAnno.

Recall from the Hypergeometric distrbition -- the analogy is that we have a bag of marbles with two different colors, say one blue and another red. A sample from such distribution consists of marbles drawn from this bag, and a probability measure for this distribution is the number of blue marbles in the sample.

Could someone explain what is the correspondence between the hypergeometric test and the significance overlap between two different set of genomic intervals?

# Answer
There is no such thing as a hypergeometric test, at least in statistical textbooks. It's a fisher test based on hypergeometric distribution.

If it is chip-seq for the same target, i.e biological replicates, significance of overlap is not quite meaningful. You get more information by for example checking the correlation of the coverage between your replicates using deeptools or anything else.

If you are comparing between two different targets, for example, chip for two different histone marks, you can use the fisher test or chi-square this way.

1. Partition the genome in X bp bins, these are all your marbles
2. Identity those that overlap with peaks from target 1 and those that don't
3. Do 2. for target2
4. construct a contingency table using 2 and 3.

So you should end up with something like this (purely made up):

```
                      no_overlap_target2 overlap_target2
  no_overlap_target1               9799             112
  overlap_target1                    88               1

```

And if there is overlap, the bottom diagonal with be higher.. less on the off-diagonal, and you can test it with a fisher.test.
> 6 votes

# Answer
Just to make sure, I took a quick look at chippeakanno, and it uses a hypergeometric test for gene enrichment, not necessarily overlap of two chip-seq results.

As you mentioned, you can imagine in your set of genes, you have those that overlap with your peaks, and those that do not. Hypergeometric distribution calculates the probability of getting a certain number of genes that overlap with your peaks (which is the same as sampling some balls and calculate the p-value of some balls with one color). The hypergeometric test calculates a p-value to see whether this number of overlapping genes is more than expected.

I also read on Wikipedia that the P-value of hypergeometric test times 2 is exactly the P-value for Fisher's exact test. Fisher's exact test tests whether your overlap is more than (or less than) expectation. Hypergeometric tests just feel like a one-tailed version of Fisher's exact tests.
> 1 votes

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


# Question
Title: Is the number of contigs and scaffolds an artifact of genome assembly and rather arbitrary, or is it canonical for each species?
E.g. bowhead whale genome has X contigs and Y scaffolds, but couldn't these easily be different numbers?

# Answer
> 1 votes
This is purely arbitrary and an artifact of how well the assembly turned out. Better quality assemblies with have fewer and longer contigs.

---
Tags: genomics
---


title: Speeding up network creation in WGCNA

text:
# Question
Title: Speeding up network creation in WGCNA
I am following the tutorials on the WGCNA website for creating a coexpression network. It is taking way too long. Is there a way to speed up execution? All I need is to create the co-expression graph. I will then export to Python and do my analysis there. In particular, the step "calculating module eigengenes" is taking a lot of time. I do not really care about the modules. All I want is the graph with the genes and the weighted edges.

# Answer
You can use the "manual network construction" approach. Calculate TOM, then use `exportNetworkToCytoscape` or `exportNetworkToVisANT` to export the network (or simply save the TOM matrix as a table and load it in Python).
> 0 votes

---
Tags: wgcna
---


# Question
Title: Get regions' information from DNA sequence data (bsgenome.hsapiens.ucsc.hg19)
I have a problem in R. I have the following dataSet (the first three rows shown) (the 5th number is the methylation level in its region):

```
chr1    15864   15866   FALSE   894 +
chr1    534241  534243  FALSE   921 -
chr1    710096  710098  FALSE   729 +

```

and then I'm trying to get the sequences from every row in the dataSet (427.000 rows in total):

```
library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
chr<-as.matrix(as.character(mydataSet[,1]))
start<-as.matrix(as.integer(as.character(mydataSet[,2]))-50)
end<-as.matrix(as.integer(as.character(mydataSet[,3]))+50)
Sequences_for_50_plus_bases<-getSeq(genome,chr,start=start,end=end)

```

I'd like to find a new fuction where I can have this dataSet, and find for its rows (which represent sequences in the human genome) if they refer to a promoter, or an enhancer or any other DNA regions' information, because I'd like to compare the differences between the regions' methylation levels.

Is there a way to do this?

# Answer
You can get specific sequences efficiently using the following packages:

```
library(GenomicFeatures)
library(Biostrings)
library(tidyverse)
library(data.table)
library(dtplyr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

Methylation_data <- tribble(
  ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
  "chr1",15864,15866,FALSE,894,"+",
  "chr1",534241,534243,FALSE,921,"-",
  "chr1",710096,710098,FALSE,729,"+") # Arrange your dataset

Methylation_GRanges <- GRanges(  # Construct a GRanges objects
  seqnames = Methylation_data %>% pull(seqnames), # First pull in the chromosome locations
  IRanges(start = (Methylation_data %>% pull(start))-50, # then the start position minus 50 nt
          end = (Methylation_data %>% pull(end))+50), # the end position plus 50 nt
  strand = Methylation_data %>% pull(strand) # the strand your sequence is located on
)
> getSeq(Hsapiens,Methylation_GRanges)
A DNAStringSet instance of length 3
width seq
[1]   103 GCTGCTGCTTCTCCAGCTTTCGCTCCTTCATGCTGCGCAGC...CTTGGCGGATGGACTCTAGCAGAGTGGCCAGCCACCGGAG
[2]   103 ACAGTTCCAATGTAATCAGAGAGAACATCACACACACACCA...ACCTGAGCAGCACTCTGCAAAGCTGTCAAGGCGGTGAAGC
[3]   103 GATAAAATAAAGCTTAGATTGGAAAAAATATTTAAGATTCT...GGAAGCTGAGTAATTGTATGTTCAAATACTTGCAAAACAT

```

To write the sequences to a file, use `writeXStringSet()`

For the second part of what you need, `GenomicFeatures` has a number of functions that obtain the annotated sequences of the whole genome by its type:

```
transcriptsBy(x, by=c("gene", "exon", "cds"), ...)
exonsBy(x, by=c("tx", "gene"), ...)
cdsBy(x, by=c("tx", "gene"), ...)
intronsByTranscript(x, ...)
fiveUTRsByTranscript(x, ...)
threeUTRsByTranscript(x, ...)

```

And this is how you use one of them to find whether your sequences overlap a particular type of DNA:

```
Find_overlap <- function(Ref_Data,
                         Overlap_Message){

  Ref_Data <- as.data.table(Ref_Data) # Convert to data.table for efficiency

  Methylation_data <- tribble(
    ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
    "chr1",15864,15866,FALSE,894,"+",
    "chr1",534241,534243,FALSE,921,"-",
    "chr1",710096,710098,FALSE,729,"+") %>%
    mutate(Overlap = "-",
           start = start-50,
           end = end+50) %>%
    as.data.table()

  Methylation_data[Ref_Data,
                   Overlap := Overlap_Message,
                   on = .(start >= start,
                          end <=  end)] # wherever the range overlaps, assign an identifier to it
}

Methylated_in_5UTR <- fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene) %>%
  Find_overlap(., "5' UTR") %>%
  print

```

seqnames  start    end Logical Methylation Level strand Overlap
1:     chr1  15814  15916   FALSE               894      +       -
2:     chr1 534191 534293   FALSE               921      -       -
3:     chr1 710046 710148   FALSE               729      +  5' UTR
```

> 1 votes

---
Tags: r, bioconductor, phylogenetics, methylation, hg19
---


# Question
Title: Weighted Network Analysis Book by Steve Horvath, taking the power of a matrix
In Horvath's *Weighted network analysis*, the author mentions that M^b, where b is a real number and M is a matrix, denotes the element-wise power. This seems odd to me. The way I was taught, M^b would mean multiplying the matrix M, b times. Should element-wise power not be M^(ob) where where o denotes the Hadamard power?

# Answer
You CAN calculate M^b for b being float number, complex number or even another matrix. As well as you can calculate exp(M), log(M), sin(M) whatever... That is very standard for math people.

Let me give two ways how you can think of it.

1) M^b = exp(b log(M)), so you need to define what is exp and log for matrices. That can be done using power series: exp(M) = 1 + M + M^2/2! + M^3/3! + ... log(1+M) = M - M^2/2 + M^3/3 - ...

2) More general way to calculate any function of matrix f(M) , looks like that: if your matrix is diagonal, then it is true that: f(diag(d\_i)) = diag( f(d\_i) ) . So the question how to reduce any matrix to diagonal one. That can be done for almost all matrices by conjugation M = C D C^{-1}. Then f(M) = C f(D) C^{-1}.
> 2 votes

---
Tags: wgcna
---


# Question
Title: What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)?
What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)?

Barycenter scores are calculated as 1 / (total distance from vertex v to all other vertices) in a strongly connected network. source

Closeness Centrality (Freeman) are calculated total distance from a node v to all the other nodes in a network. source

It looks like the same centrality .

Ps: I want to analyze closeness for my protein–protein interaction network.

# Answer
> 1 votes
First, I'll observe that Freeman centrality as you define it is the reciprocal of the Barycenter centrality as you define. It looks like people often use the reciprocal instead of the original definition, so in that sense they do look the same.

I found the following comparison here (though it doesn't name it as Freeman):

> There are 2 types of distance centrality scores, Closeness Centrality and Barycenter Centrality. Barycenter Centrality for vertex v defined as:
>
> 1 / (total distance from v to all other vertices)
>
> Closeness scores are calculated using the formula 1 / (average distance from vertex v to all other vertices) and Barycenter scores are calculated as 1 / (total distance from vertex v to all other vertices).

So it appears that for this package it is a difference of the reciprocal of the total vs. reciprocal of the average.

Notably, wikipedia differs here, in that it matches your source:

> Freeman's closeness centrality, the total geodesic distance from a given vertex to all other vertices, is the best known example.citation

Yet another source notes the existence of a "normalized" version of Freeman's closeness centrality, which normalizes centrality of a node by the degree of the network. It also comments on whether it is the reciprocal or not:

> Closeness (Freeman '79) The graph-theoretic distance of a given node to all other nodes. The sum of the rows/columns of the geodesic distance matrix of a graph.
>
> Simple closeness is an inverse measure of centrality: the larger the numbers, the more distant an actor is, and the less central. Should really be called "farness".
>
> Normalized version divides the minimum "farness" possible (N-1) by "farness" to simultaneously make the range 0 to 1 and invert the measure so that larger values correspond to greater centrality (truly "closeness")

So I guess that it might depend on what tools you are using. Most of the time, it might be the same. Other times, it might be normalized in different ways with either averages or for 0-1 ranges.

Closeness centrality is probably the more standard term, so it is unlikely that anyone would criticize it. But it appears that there is sometimes some variation in what people mean by that term and whether they are using a normalized or a raw closeness centrality.

In my opinion the normalization leading to a range from 0-1 is nice because it makes the measure at least somewhat generalizable across networks. But it's up to what you need to do.

---
Tags: networks, protein-protein-interaction
---


# Question
Title: How can I convert codon coordinates to genomic position?
I am looking for a given mutation in IGV, which accepts coordinates in the form of

```
chr<X>:<Y>

```

Where X is the chromosome number and Y is the base location in that chromosome.

Eg. `chr6:26091179` which has the following attributes:

```
Gene: HFE|LOC108783645
Marker: rs1799945
Position: chr6:26091179

```

Additionally, an interval can be written like

```
chr<X>:<Y>-<Z>

```

So basically from base located at Y to Z.

**My questions are:**

**1)** How do I convert `HFE|LOC108783645` to `chr6:26091179` (and vice-versa)?

**2)** How do I convert `rs1799945` to `chr6:26091179` (and vice-versa)?

**3)** On a different example, how do I convert something like `codon 109 of exon 3 of beta globin gene (GTG->TTG)` to `chr<X>:<Y>` (or `chr<X>:<Y>-<Z>`)?

Why are there so many ways to specify positions in genes? For me the most clear way is `chr<X>:<Y>`, all the others leave me guessing...

# Answer
> 1 votes
1. `HFE|LOC108783645` is https://www.ncbi.nlm.nih.gov/gene/108783645. You can use biomaRt to go from the `LOCxxxx` ID to location. Searching by `HFE` will give you a different location though, so you might have to account for such eventualities. It is a long and often manually intensive process.
2. You should be able to use entrez-direct/eutils to go from rsID to position
3. Unless there's some sort of natural language processing API, `codon 109 of exon 3 of beta globin gene` cannot be parsed by existing APIs (Unless it's always `codon X of exon Y of Z`, in which case you could use `grep` too get the components. Even so, going from exon + codon + codon base change to exact base location will require a lot of wrangling.

---
Tags: gene, genome, phylogenetics, igv
---


title: phangorn fasta file read error

text:
# Question
Title: phangorn fasta file read error
Im trying to read a multifasta sequence file into phangorn I get the following error

```
tipseqfile <- ("abhi_seq/clean_dup_rem.fasta")
tipseq <- read.phyDat(tipseqfile,format="fasta")

```

my file

> Error in read.phyDat(tipseqfile, format = "fasta") : sequences have different length

Now is it necessary to have same length of each fasta sequence? or Im doing something wrong

# Answer
> 3 votes
You have to perform an alignment first. Phangorn use phylogenetic estimation and this requires a precise alignment position to track homologous nucleotide mutation. Maximum likelihood theory will assess the probability of a given tree for each homologous nucleotide site, if they are not alignment it can't do it and will violate the assumptions of phylogeny.

If the sequences are different lengths, which it means both nucleotide and indel ('-') sites, it means there hasn't been or there is an error with the alignment. I use `muscle` these days having previously used `clustalo`. Its really easy to use and easily configures on a linux server and its super quick. If you get stuck let me know.

I had a look at your alignment and can verify it hasn't been aligned.

---
Tags: r, phylogenetics, fasta, phylogeny, phangorn
---


# Question
Title: Inflated p-values in quantitative trait analysis
I am performing a quantitative trait association between the expression of one gene and ~400,000 methylation values. First, both variables are rank inverse normal transformed, adjusted for confounders (technical and for family structure) in a mixed-effects model and we use the residuals for the association. Then we do a linear regression between the residuals of expression and each methylation site, adding other four confounders in the model (age, BMI, etc.). The sample size is of ~500. For the association we are using MatrixEQTL.

When we check the p-values, they are clearly inflated. Is it something that I have to be concerned about? Did I probably miss something in the analysis? Thanks in advance.

# Answer
Permutation as suggested by @StupidWolf's comment is essential to understand what's going on. If permutation makes this pattern go away, then you have a problem with your model specification, there's something uncorrected.

If your data are weird, well, that's just how they are. But this argues to me that something else is going on confounding your associations. Cryptic population structure in the data would lead to this kind of pattern, e.g. here. That is the first thing that most reviewers would say, I'd guess.

Some random thoughts:

1. What does the distribution of your methylation data look like after normalization?
2. Your observed distribution is basically never on the expected line, which argues strongly for confounding.
3. How are you adjusting for non-independence of methylation types? **Family-wise correction is likely not specific enough. You may have to estimate an identity-by-state SNP similarity matrix for your population and fit that as part of your model,** or at least fit 10-20 principal components to account for population structure.
4. How are you representing your methylation status? If it's binary you might look at this paper.
> 4 votes

---
Tags: statistics, gene-expression, methylation, linear-regression, qtls
---


# Question
Title: Receptor-Ligand R-package for visualization
I was wondering if there was any R-package for receptor-ligand interaction similar to SingleCellSignal (https://github.com/SCA-IRCM/SingleCellSignalR\_v1/blob/master/SingleCellSignalR/vignettes/UsersGuide.Rmd)

+++ if there are vignettes and/or the graphs can be customized.

+++ if it can be linked to Seurat

Note: I have tried iTALK, but I keep getting an error that seems to be a common problem and have left a query.

Thanks

# Answer
> 1 votes
Two nice tools for receptor-ligand interaction analysis are:

* CellPhoneDB by the Teichmann lab, although this is not an R package...
* NicheNet by the Saeys group

Note that these tools are mainly meant for analysis of receptor-ligand interactions, although I think NicheNet also has some visualization options. CellPhoneDB needs to be run from a shell, but the results can easily be imported into R for visualization. You might check the papers on these 2 methods for some ideas of visualization. Combining Seurat's visualization options and ggplot2 should get you already pretty far.

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