# Question
Title: pvalue calculation of protein-protein network with permutation test
*Note: this question has also been asked on Biostars*

I have performed a network recreation analysis based on the interactions of proteins from String db.

I wanted to know if the interactions have significance by reshuffling the proteins with permutation against the background network. Since, its first time for me to perform permutation test on PPI, Can someone guide me to a tool or a tutorial to perform permutation test on newly created PPI network.

I tried looking at a youtube tutorial and read some articles like this one, but they were not very useful for my purpose.

# Answer
It is probably necessary to correct for the degree distribution of the network, I am not sure that a tool such as the Mantel test will quite do this (as suggested by user Michael G.). In naturally occurring networks degree distribution is a strong confounder.

Then again I am not totally sure what shape your data is in or what the question is. You mention a background network, but I am not sure what this is or how you want to use it. Possibly your network reconstruction is a subset of your background network? In which case you just define the rule for how your subsample is drawn and then you compare it to permuted networks.

However, here is some fairly crude but serviceable R code that implements a general purpose network rewiring test, using igraph. I am not using a background network as again I'm not sure what it is.

(Python also has igraph so should be similar `rewire()` function)

```
compute_measure = function(network) {
  # i don't know what you are measuring, so you have to write this function.
  # this function takes a graph object as input
  # this function returns a single numeric value (whatever you're measuring) as output
}

# "net" is a network of the graph type from igraph

measure = compute_measure(net)

num_perms = 1000

perm_measures = c()

for (i in 1:num_perms) {
  # permute the network by rewiring, rewire 10x as many times as edges
  # to ensure saturation
  perm_net = rewire(net, with = keeping_degseq(niter = gsize(net) * 10))
  perm_measures[i] = compute_measure(perm_net)
}

# one-sided monte carlo p value is proportion of rewirings (plus one) greater than your measure
p_value = (1 + length(which(perm_measures >= measure))) / num_perms

```

If on the other hand you are just trying to come up with a single permutation of the edges of your network to compare against, you can run `rewire()` only once to get that permuted network. It depends on what question you are trying to answer, and what null hypothesis you're trying to get a p value for.
> 1 votes

# Answer
**Two answers:** What you have asked for is a straight Mantel's test. You could read up on the test, it is traditionally used to assess genetic vs geographic distance but could be applied. This is easy to do in R:

```
mantel.test(m1, m2, nperm = 999, graph = FALSE,
            alternative = "two.sided",  ...)

```

You will need to import your data (recommend csv) and parse to a data.frame both for your observed (m1) and expected (m2). Its easy.

The documentation is here. My advice is to think carefully about the contruction of the matrix: its an everything versus everything comparison.

**Alternative approach** My personal choice would be to use a covariance/correlation matrix. That can be done in R too, but you need to understand what covariance is and how it relates to your problem. If you don't really get covariance, I'd leave it. Hence I'll avoid putting code down - but is part of ggplot2.

Please note, I don't supply detailed codes (I'm not your supervisor), unless its Perl, which this isn't
> 0 votes

---
Tags: python, statistics, protein-protein-interaction
---


# Question
Title: Save output to a specific folder and/or with a specific prefix in Cancer Genomics Cloud
I would like to run the CNVkit workflow on Cancer Genomics Cloud.

However, the two workflows that I'd like to use:

don't have parameters that can be set to define in which folder to save the results or which prefix to add to the results. This means that if you're generating multiple references or calling CNV for multiple tumor samples you will end up with results having the same name for each run (tumor).

I tried to modify the CWL workflow, online and through Rabix composer, but I was unable to make it work. I would like to use the published workflow to not code everything from the ground up.

Do you know how to update/use the current CWL workflow in order to avoid having the same filename/folder for all the results?

# Answer
If the output is File, you can add this step https://github.com/Barski-lab/workflows/blob/master/tools/rename.cwl to rename the output file.
> 1 votes

---
Tags: cnv
---


# Question
Title: hgsql not found when running a script despite being accessible on command-line
I am working on macOS. I am following the steps described here to lift over an annotation from the reference genome of *C. elegans* to another assembly of *C. elegans*.

I have run this code:

```
time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=`pwd` \
      -ooc=`pwd`/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \
  -bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \
  -target2Bit=`pwd`/${target}.2bit -targetSizes=`pwd`/${target}.chrom.sizes \
  -query2Bit=../${query}/${query}.2bit \
  -querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1

```

However, the `do.log` file contains:

```
zsh:1: command not found: hgsql
HgStepManager: executing from step 'align' through step 'cleanup'.
HgStepManager: executing step 'align' Sun Nov 11 15:38:41 2018.
Using localhost, /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit and ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit
# chmod a+x /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/job.csh
# chmod a+x /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh
# ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' localhost nice /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh
cd /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat
rm -rf tParts
/Users/username/Documents/projects/main/output/software/liftover/data/scripts/partitionSequence.pl 58018742 0 /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.chrom.sizes 2000 -lstDir=tParts
rm -rf qParts
/Users/username/Documents/projects/main/output/software/liftover/data/scripts/partitionSequence.pl 10000000 0 ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.chrom.sizes 1000 -lstDir=qParts
Can't open seq size file ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.chrom.sizes: No such file or directory
Command failed:
ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' localhost nice /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh

```

I think the biggest problem is that the `hgsql` command is not found.

If I run `hgsql` on the command line, I find it even though I get the following error:

```
dyld: Library not loaded: /usr/local/opt/mysql/lib/libmysqlclient.20.dylib
  Referenced from: /usr/local/bin/hgsql
  Reason: image not found
[1]    34059 abort      hgsql

```

Location of the hgsql installation:

```
which hgsql
/Users/username/Documents/projects/main/output/software/liftover/data/bin/hgsql

```

Can anyone help me with using hgsql, please?

# Answer
It is possible that the shell being used in the process is different from your standard command line shell. I am not familiar with the tool, and it's hard to say more without knowing your OSX version, but it looks like it is using zsh. I think that zsh only fairly recently became the standard OSX shell. I don't know the tool or the shell very well, but it looks like you might be able to explicitly set up hgsql in your PATH in your .zshrc.

It's also possible that the liftover script was not built to run via zsh. I notice that it is trying to run csh scripts as part of the workflow. Then again the localhost may insulate it.

So I think that my first step would be to sort out which shell I am using, whether it works with the liftover pipeline, and how it's running on the localhost. I notice that the ssh call running the csh script on the localhost is the step that actually fails. What shell is executing the csh script? Does that shell know about hgsql? It could be that the way that the error bubbles up it is hiding which shell it is actually using.

I'd recommend trying that ssh call in isolation from the rest of the pipeline, and diagnosing the error that way.

Hope that helps.
> 1 votes

---
Tags: sequence-annotation, software-installation, errors, liftover
---


title: design formula question

text:
# Question
Title: design formula question
I am not sure I am building the proper design formula for the question I want to answer

I have the following samples with three factors; clone, the structure and the condition.

```
clone structure diabetic
1 07        2D      Dia
2 21        2D     Ctrl
3 23        2D      Dia
4 32        2D     Ctrl
5 34        2D      Dia
6 43        2D     Ctrl
7 07        3D      Dia
8 21        3D     Ctrl
9 23        3D      Dia
10 32        3D     Ctrl
11 34        3D      Dia
12 43        3D     Ctrl

```

I want to pull the differentially expressed genes for structure (2D vs 3D) blocking by clone, and the one differentially expressed by condition (regardless of the structure variable).

I am not sure how to build the design in one call, therefore I have produced the first object to pull 2D vs 3D blocking by clone using the following.

```
des_structure <- DESeqDataSetFromHTSeqCount(sampleTable = samples,
                                       directory = "../data/htseq_geneCounts/",
                                       design= ~ clone+structure)

```

for the condition controlled by structure I have used

```
des_condidtion <- ddsHTSeq_diabetic <- DESeqDataSetFromHTSeqCount(sampleTable = samples,
                                                 directory = "../data/htseq_geneCounts/",
                                                 design= ~ structure+diabetic)

```

I am not sure whether this is the best way to build the designs.

# Answer
> 6 votes
You can't (easily) use a single design, since clone is nested within diabetic status, so you were correct using two separate designs. Your designs are correct.

If you really wanted to use a single design for some reason, it'd be `~clone+structure` and for the latter comparison you'd end up using a rather complicated contrast (something like `c(rep(c(1,-1), 3), 0)`). I really can't recommend that unless you have a very compelling reason, since it's vastly easier to think in terms of factors for your experiment and you're not hurting anything by fitting twice.

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


# Question
Title: significant SNPs to annotated candidate genes
I have performed GWAS and got some significant SNPs. How can I get candidates genes in the flanking regions of those SNPs and also their functional annotations (eg. KEGG).

# Answer
You can use BEDOPS' closest-features tool. You only need to have the SNPs coordinates in .bed format and a .bed file with gene locations. BEDOPS will look for the closest genes upstream and downstream of your SNPs. The command will look something like this:

```
closest-features yourSNPs.bed yourAnnotation.bed

```

Then, you can obtain your functional annotations through different tools such as biomart.
> 2 votes

---
Tags: snp, gwas
---


# Question
Title: How do I get the adjacency matrix from gene similarity network?
I am exploring application of Graph convolution network on gene similarity network. I am trying to explore genemania.org and https://string-db.org/

Please let me know where can I find the data for my research in gene interaction network, which allows me to download the data in the form of an adjacency matrix.

# Answer
The STRING download page is here. I am not sure that they provide any data in the form of an adjacency matrix, as the number of nodes (rows/columns) would be quite large, and for this reason they are likely to represent the data as lists of links, for compactness.

Their data will likely have to be cleaned for any purpose you are interested in. You will have to make some decisions on how to filter the data. For example, as suggested in the comments, you may want to filter on the quality of the evidence supporting interactions, or based on species of the protein, or other things.

It's possible that you may have better luck going with another database that has done some of these steps already. For example, there are relatively compact human interactome datasets here, or yeast interactome datasets here. In this case you will still need to convert from a list of links to an adjacency matrix, however this is not very difficult if you can use a scripting language.

Notably, you may have to figure out how to distinguish missing data from tested non-interactions in these cases (or decide that missing data is ok). In most cases people will not have exhaustively tested every cell of the all-by-all interactome adjacency matrix.

Hope that helps.
> 1 votes

---
Tags: gene, networks
---


# Question
Title: Single sample in group: normal pipeline or Kal's Z test
As stated, which one is better for differential expression analysis? When I say normal pipeline I mean limma-voom, edgeR and DESeq2 pipeline for standard analysis. Kal's z test is mentioned in this paper and I know some pre-made tools like the one on QIAgen uses this as well for single vs single comparison

The situation I have is that I have RNA-Seq data for some patients with rare translocations in leukaemia. There are two of these translocation with almost 10 or even more patients so I have no problem with them. And then there are these translocation which mostly I have only 1, some with 2 or 3 patients, for 7 different translocations.

Our hypothesis is that the two larger group can certainly form their own subtypes, and the others may fall into one of the three cases: (1) belongs to one of the two groups despite of the difference in translocation; (2) belongs to a separate group on their own i.e. complete outlier; (3) belongs to a separate group with other members in the "others" i.e. forming a third subtype. And this can be tested with clustering.

For that, I need to do supervised feature selection to filter out useless genes. I can do that with those with multiple samples, but not those isolated cases. So I thought I can get some DGE analysis and use a certain adj.p and fold change as cut off, or adj.p as ranking matrix, and select those genes for clustering.

So my main questions are: (1) when using limma/edgeR/DESeq2, should the model design include all the translocation labels in one go or should I just add the one contrast I want to make and build a new model everytime? (2) Is Kal's z test useful for one vs all DGE analysis and is it better? (3) Are there other outlier analysis working better under this situation?

Many thanks!

# Answer
1. You'll have a bit more power if you use all of the samples at once.
2. Nothing is useful for this, don't waste your time with it.
3. You're not looking for outliers, you're looking to see how samples cluster. You can just use correlation for that.
> 1 votes

---
Tags: differential-expression, deseq2, edger, limma
---


# Question
Title: How are residuals used as new phenotype?
I'm a novice trying to wrap my head around with fancy jargons. If I search residual, I understand that it represents the difference between the observed value of the dependent variable (y) and the predicted value (ŷ).

I'm unable to follow then, how are these residuals used as new phenotype?

I found similar question in the forum: Adjusting phenotypes by regressing out covariates

The origin of my question is predixcan pipeline, where residuals are used as the new gene expression values after regressing input values on covariates.

# Answer
> 1 votes
The residual is the level of error in a regression model, the lower the residual the better the model. Residuals cannot equate to a phenotype, it is actually the opposite if your regression model is trying to investigate a given genotype/phenotype(s) model. The residual in this case is the amount of infomation which the genotype/phentype does not account for, i.e. there are other factors resulting in your observations which you have not captured in your model.

What I assume is that the regression model you are referring is investigating whether the residual is lowered by examining different gene combinations. The objective of a regression model is to lower the residual and I suspect that it is comparing the new residual against the old one. Therefore it is not being used as an input into a regression model but an assessment of its robustness.

---
Tags: statistics, linear-regression
---


# Question
Title: Correlate DEGs from DESeq2, EdgeR and Limma results
I have a lists of DEGs identified by DESeq2, EdgeR and Limma. I would like to correlate the the gene rankings in the lists to decide on a package to use in downstream analysis. I am havig a few problems with this, would anyone be able to help?

Thanks

# Answer
I would use the p-value/FDR which each method returns to rank the gene list for that method in order from 'most likely to be DE' to least likely. You would then have three ranked lists of genes - you would put the actual rankings in for each gene - 1 for the top ranked, 2 for the next, etc. You could then carry out, say Spearman's rank correlation on the lists of rankings to assess agreement between the methods.

You could also decide to keep all 3 methods, and accept genes as DE if, for example, at least 2 out of 3 methods give an acceptable p-value/FDR.

From my experience, correlation should be reasonably high at the top of the list, where the methods tend to agree more on the rankings, and when the correlation gets lower, further down the list, would be where I would start to question whether the genes are really differentially expressed.
> 1 votes

---
Tags: rna-seq, deseq2, edger, limma, correlation
---


# Question
Title: Appropriate tool or algorithm for sloppy alignment of degenerate bases
I have an optimization problem where I have a degenerate nucleotide sequence I want to align against subsets of a reference genome (exons, specifically, to make the problem more tractable).

The sequence itself contains subregions of variable length.

The sequence looks like this, containing degenerate bases:

```
           L           C           R
5' - SSS......SSS|NNN.....NNN|WWW......WWW - 3'

```

Per IUPAC notation, `S` denotes `C` or `G`, `N` denotes any base, and `W` denotes `A` or `T`.

I am marking each block `L`, `C`, and `R` to broadly denote "left", "center", and "right" regions, oriented from 5' to 3' ("forward"-strand oriented).

The widths of the `L` and `R` blocks can be from 20 to 25 bases long. The width of the `C` block can be 14 to 16 bases long.

One programmatic way I could think to do this would be to brute-force a solution, generating all combinations of sequences of varied lengths, `blastn`-aligning each of such combinations to the sequence for the exon, and then looking for the best alignments among such sequences.

The output would be the scores for regions from sequences that align, to be able to say, from point X to point Y of the genome, a sequence here aligns better than other candidate sequences.

To avoid reinventing tools, and avoid brute-force methods, I'd like to know if there is a tool or even a wrapper around tools, which does this, in case it is already written, or another approach that is more elegant/efficient.

# Answer
> 1 votes
I believe mafft can align degenerate sequences. From the documentation for version 7 at: https://mafft.cbrc.jp/alignment/software/anysymbol.html

```
Acceptable nucleotide symbols in the default mode:
a, c, g, t, u
Ambiguous nucleotides (r, y, w, s, k, m, d, v, h, b; IUPAC-IUB codes) can be used and are scored as:
score(r,a) = ( score(a,a) + score(g,a) ) / 2
score(r,g) = ( score(a,g) + score(g,g) ) / 2
score(r,t) = ( score(a,t) + score(g,t) ) / 2
score(r,c) = ( score(a,c) + score(g,c) ) / 2
...
score(r,y) = ( score(r,t) + score(r,c) ) / 2
...
score(r,r) = ( score(a,a) + score(g,g) ) / 2
score(y,y) = ( score(c,c) + score(t,t) ) / 2

```

MAFFT is designed for multiple sequence alignment, so will attempt to find a local alignment of homologous segments (using FFT), and then attempt to make a global alignment between the input sequences.

---
Tags: sequence-alignment, blast, sequence, blastn
---


# Question
Title: DNA bioinformatic search and methylation analysis by qMSP
Im quoting the methods which was used

> We performed a bioinformatics analysis to identify potential CpG islands in the genomic sequences of RANKL and OPG genes. The USCS Human Genome Browser public database (genome.ucsc.edu), the Primer Express software (available from Applied Biosystems), and CpG Island Explorer software (cpgie.sourceforge.net) were used to explore the nucleo- tide sequences surrounding de transcription start site (TSS) of RANKL isoforms I and II (NG\_008990.1 Ref Seq Gene), and OPG (NG\_012202.1 Ref Seq Gene). A (C + G)/total bases ratio
>
> > 0.5 and CpG observed/CpG expected ratio \>0.6 were used as criteria.

Im looking for the region which they have used

This is the paper

```
"Role of DNA methylation in the regulation of the
RANKL-OPG system in human bone"

```

I tried using ucsc to check the region but didn't get much. The idea is to find the region to design the primer for the same

How can i get the region from ucsc ,the CpG island around the TSS of the above mentioned genes which are RANKL and OPG

Any suggestion or help would be highly appreciated

# Answer
> 1 votes
It's likely that they accomplished this by using the UCSC Table Browser to:

1. Find the loci of the two isoforms you mention
2. Finding CpG islands (Group = Regulation, Track = CpG Islands) near those transcription start sites
3. Design primers using Applied Biosystem's software

Alternatively, they list the exact primers used in Table 1.

They don't detail exactly how they did it, so if you have specific questions it might be worth it to email the corresponding author

---
Tags: methylation
---


# Question
Title: How to plot character state changes from a presence-absence matrix on to a phylogeny
I wish to assign character state changes from a presence-absence matrix to a phylogeny. Thus I want to identify the most parsimonious hypothesis for which node and branch a given "mutation" or "change" has occurred.

I have tried assigning each character to its leaf node, and then if the leaf node's sister has the same character I reassign the character to the parent node (and work back until all nodes are assigned). I am using a dummy dataset to try to achieve this:

```
Matrix
>Dme_001
1110000000000111
>Dme_002
1110000000000011
>Cfa_001
0110000000000011
>Mms_001
0110000000000011
>Hsa_001
0110000000000010
>Ptr_002
0110000000000011
>Mmu_002
0110000000000011
>Hsa_002
0110000000000011
>Ptr_001
0110000000000011
>Mmu_001
0110000000000011

Phylogeny
((Dme_001,Dme_002),(((Cfa_001,Mms_001),((Hsa_001,Ptr_001),Mmu_001)),(Ptr_002,(Hsa_002,Mmu_002))));

```

I assign internal nodes using ete3, so my output should be:

```
BranchID    CharacterState    Change
Node_1:    0    0->1
Hsa_001:    15    1->0

```

As my code assigns character states based on their sisters if a loss is encountered it messes up the output so that:

```
BranchID    CharacterState    Change
Node_1:   0     0->1
Node_3    15    0->1
Node_5    15    0->1
Node_8    15    0->1

```

I'm coding in python 2.7 and welcome assistance.

My code:

```
from ete3 import PhyloTree
from collections import Counter
import itertools

PAM = open('PAM','r')

gene_tree = '((Dme_001,Dme_002),(((Cfa_001,Mms_001),((Hsa_001,Ptr_001),Mmu_001)),(Ptr_002,(Hsa_002,Mmu_002))));'

NodeIDs = []

tree = PhyloTree(gene_tree)
edge = 0
for node in tree.traverse():
    if not node.is_leaf():
        node.name = "Node_%d" %edge
        edge +=1
        NodeIDs.append(node.name)
    if node.is_leaf():
        NodeIDs.append(node.name)

f = open('PAM','r')

taxa = []
pap = []

for line in f:
    term = line.strip().split('\t')
    taxa.append(term[0])
    p = [p for p in term[1]]
    pap.append(p)

statesD = dict(zip(taxa, pap))

def PlotCharacterStates():

    Plots = []

    events = []

    for key, value in statesD.iteritems():
        count = -1
        for s in value:
            count+=1
            if s == CharacterState:
                a = key, count
                events.append(a)

    Round3_events = []
    while len(events) > 0:
        for rel in Relationships:
            node_store = []
            sis_store = []
            for event in events:
                if rel[0] == event[0]:
                    node_store.append(event[1])
                if rel[1] == event[0]:
                    sis_store.append(event[1])
            if (len(node_store) > 0) and (len(sis_store) > 0):
                place = rel, node_store, sis_store
                Round3_events.append(place)

        moved = []
        for placement in Round3_events:
            intercept = (set(placement[1]) & set(placement[2]))
            node_plot = (set(placement[1]) - set(placement[2]))
            sis_plot = (set(placement[2]) - set(placement[1]))
            if len(node_plot) > 0:
                for x in node_plot:
                    y = placement[0][0], x
                    Plots.append(y)
                    moved.append(y)
            if len(sis_plot) > 0:
                for x in sis_plot:
                    y = placement[0][1], x
                    Plots.append(y)
                    moved.append(y)
            if len(intercept) > 0:
                for x in intercept:
                    y = placement[0][2], x
                    y1 = placement[0][0], x
                    y2 = placement[0][1], x
                    moved.append(y1)
                    moved.append(y2)
                    events.append(y)

        for event in events:
            if event[0] == "Node_0":
                Plots.append(event)
                moved.append(event)

        events2 = (set(events) - set(moved))
        events = []
        for event in events2:
            events.append(event)

    pl = set(Plots)
    Plots = []
    for p in pl:
        Plots.append(p)

    print CharacterState, Plots

'''
assign sisters to leaves, internals
'''

e = []
round1b_e = []
round2a_e = []
placements = []
Relationships = []
Rounds = []
for node in tree.traverse():
    sisters = node.get_sisters()
    parent = node.up
    cycle1 = []
    if node.is_leaf():
        for sister in sisters:
            if sister.is_leaf():
                round1a = ["Round1a", node.name, sister.name, parent.name]
                node_names = node.name, sister.name
                Rounds.append(round1a)
                e.append(node_names)
                x = node.name, sister.name, parent.name, "leaf-leaf"
                Relationships.append(x)
            if not sister.is_leaf():
                round1b =  ["Round1b", node.name, sister.name, parent.name]
                node_names = node.name, sister.name
                Rounds.append(round1b)
                round1b_e.append(node_names)
                x = node.name, sister.name, parent.name, "node-leaf"
                Relationships.append(x)
    elif not node.is_leaf():
        if not node.is_root():
            for sister in sisters:
                if not sister.is_leaf():
                    node_names = node.name, sister.name
                    round2a_e.append(node_names)
                    x = node.name, sister.name, parent.name, "node-node"
                    Relationships.append(x)

x = []
CharacterStates = []
for key, value in statesD.iteritems():
    for value in value:
        x.append(value)

y = sorted(set(x))
for x in y:
    CharacterStates.append(x)

for CharacterState in CharacterStates:
    PlotCharacterStates()

```

# Answer
To give some formal words for what you are doing, I would say this sounds like you are trying to reconstruct ancestral states. There has been some work done on this, and a similar algorithm to yours was devised by Fitch in the 1970s: Fitch W: Toward defining the course of evolution: minimum change for a specific tree topology. Syst Zool 1971, 20: 406–416. 10.2307/2412116

The missing part, from your approach, is that once Fitch gets to the root of the tree, he than retraverses back to the leaves, resolving ambiguities along the way (if two sister nodes have different states, how do you decide what the parental state was? Fitch's idea is that once you get to the root, you can come back down and use the parental state to set the child state in ambiguous cases).

There is a nice implementation of Fitch's algorithm in the R phangorn package, which makes the reconstruction a one-liner (given a rooted tree object, and a matrix of leaf node states):

```
anc.mpr = ancestral.pars(rooted_tree, tree_samples_data_matrix.phyDat, "MPR")

```

If you don't want to use the R package itself, you could perhaps use it as pseudocode to port into python.
> 1 votes

# Answer
The OP could use the Rpy2 library to connect Python to R, but it is often better to use R directly.

The alternative is to simply use MacClade 4.0, although you will need an old Mac. This is a GUI interface that will allow you to identfy the parsimonious character state for each node in the phylogeny. The API allows you to interactively explore the phylogeny by point and click moving the branches around the tree. The tree will then be assessed in terms of "number character steps" and automatically calculate this for each topology loaded. The output is very colourful so its easy to understand a large data set.
> 0 votes

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


# Question
Title: Visually summarizing phylogenetic networks
After having worked with Bayesian phylogenetic tree inference for some time, I am now trying out a method to infer network phylogenies. To get an idea of what it infers, I'd like to find a way of summarizing the results.

The output of the inference is a NEXUS file containing a stream of networks encoded in eXtended Newick as trees with duplicate node names, like this:

```
(((((((Kiraman:0.028228,Kui:0.028228)S1:0.0797301)#H4:0.35229)#H2:23.4168)#H1:15.8027,(Kamang:29.3316,#H1:5.45455[&&NHX:gamma=0.36825381342449554])S11:10.3482)S12:1.27331,((Klon:1.16727,((#H4:0.524833[&&NHX:gamma=0.33697036713230577],Kafoa:0.632791)S6:0.0615575)#H3:0.472921[&&NHX:gamma=0.6642547809404831])S2:35.6082,(#H2:22.2107[&&NHX:gamma=0.7544913491697585],((A-Petleng:0.421975,(A-Takalelang:0.168615,(A-Ulaga:0.145472,A-Atimelang:0.145472)S4:0.0231422)S5:0.253361)S7:0.902144,#H3:0.629771)S3:21.3468)S8:14.1045)S9:4.17761)S10:54.5301):0;

```

There is no formal restriction on the shape of the network, and hybridization edges have a weight 0\<γ\<1. How do visually or formally summarize a stack of several thousands of these networks, and in particular the backbone tree structure plus a glimpse of the hybridization events they display?

I have played with some tree summarizing approaches and tried to first infer a summary tree and then glue some cross-edges in there, but the way I came up with did not look reasonable.

I am now considering to calculate distances between leaves and using network visualization tools like SplitsTree or NeighborNet to get a visual representation of that, but how shoud I represent the uncertainty of the distances in even a single network (not to mention thousands of them) in the input for these algorithms?

Are there any good approaches to summarize phylogenetic networks of the shape given above? Or to consistently simplify them to a shape that can be easily summarized?

# Answer
> 1 votes
Firstly you treat this as two separate independent calculations.

1. The traditional approach is simply to produce a consensus phylogeny using your Bayesian method and assess the equivocal bootstraps/Bayesian probability.
2. Feed your data set (not the trees) into NeighborNet within SplitsTrees. There are statistics that classify the level of departure from a bifurcating phylogeny.

Examine both outputs for orthogonal support for incongruence.

**Note** Phylogenetic network statistics are a relatively soft area of phylogeny (they are not used as hard evidence).

**Networks from Bayesian trees**

The final approach (which is what you are asking about), is to feed your Bayesian trees into Dendroscope, this is a GUI that represents the uncertainty between different bifurcating phylogenies as networks and you can use it to understand the incongruence between different trees within your Bayesian MCMC, because the phylogenetic networks directly corresponds to the level of incongruence.

Dendroscope is a nice Java API and is very easy to use.

---
Tags: visualization, phylogenetics, networks, phylogeny
---


title: Merging two .txt files

text:
# Question
Title: Merging two .txt files
I have two .txt files and let's call it '**File A**' and '**File B**.Here are examples of what file A and file B looks like:

**FILE A**

```
#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno A
2   2000    2000    T   G   Pheno F, Pheno G
3   3000    3000    C   G   Pheno L

```

**FILE B**

```
#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno B
2   2000    2000    T   G   Pheno H, Pheno I
3   3000    3000    C   G   Pheno M
3   3000    3000    C   G   Pheno N

```

Here's how I want the output to look like:

```
#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno A, Pheno B
2   2000    2000    T   G   Pheno F, Pheno G, Pheno H, Pheno I
3   3000    3000    C   G   Pheno L, Pheno M, Pheno N

```

Both are tab-separated .txt files and I am trying to merge them together for SNP annotation and would love to hear suggestions on how to proceed this.

# Answer
Whereas wouldn't be as fast as `awk`, here is a rather dirty R solution:

```
library(dplyr)
library(tidyr)
library(magrittr)

# File 1

# the "fill" parameter is used since not all the rows have the same number of cols
file_1 <- read.table("file_1.txt", skip = 1, fill = TRUE)

# paste() function with argument collapse is applied to the "phenotype"
# columns of each row (these cols start at position 6 and goes until
# the number of cols: dim(x)[2] gives the number of cols of a table
phenotypes <- apply(file_1[, 6:dim(file_1)[2]], 1, function(x) paste(x, collapse = ""))

# Cols referring to the genomic position, ref and alteration (cols 1:5)
# are combined with the newly created "phenotypes" column
file_1 <- data.frame(file_1[, 1:5], phenotypes)

# colnames are added, could be avoided if dealt with at file reading step
names(file_1) <- c("Chr", "Start", "End", "Ref", "Alt", "Phenotypes")

# File 2

file_2 <- read.table("file_2.txt", skip = 1, fill = TRUE)

phenotypes <-apply(file_2[, 6:dim(file_2)[2]], 1, function(x) paste(x, collapse = ""))

file_2 <- data.frame(file_2[, 1:5], phenotypes)

names(file_2) <- c("Chr", "Start", "End", "Ref", "Alt", "Phenotypes")

# The two tables generated as above are merged using the appropriate
# columns as "anchors" in the "by" argument of full_join()
# The resulting Phenotypes.x and Phenotypes.y cols are concatenated
merged_tables <- full_join(file_1,
                           file_2,
                           by = c("Chr", "Start", "End", "Ref", "Alt"))
merged_tables$Phenotypes <- paste(merged_tables$Phenotypes.x,
                                  merged_tables$Phenotypes.y,
                                  sep = ",")
merged_tables$Phenotypes.x <- NULL
merged_tables$Phenotypes.y <- NULL

# Duplicated rows (in terms of genomic pos, ref and alteration) in any
# of the files result in duplicated rows in the merged table and are
# "handled" by aggregating/summarizing the info at the phenotype cols
# by using paste()
# strsplit() is used to "split" the "phenotypes" on "," in order to
# remove repeating phenotypes resulting from the full_join() call

merged_tables %<>% group_by(Chr, Start, End, Ref, Alt) %>%
      summarize(Phenotypes = paste(Phenotypes, collapse = ",")) %>%
      ungroup() %>%
      mutate(Phenotypes = strsplit(Phenotypes, ","))

merged_tables$Phenotypes <- lapply(merged_tables$Phenotypes, function(x) paste(unique(unlist(x)), collapse = ","))
merged_tables$Phenotypes <- unlist(merged_tables$Phenotypes)

# A tibble: 3 x 6
    Chr Start   End Ref   Alt   Phenotypes
  <int> <int> <int> <fct> <fct> <chr>
1     1  1000  1000 A     T     PhenoA,PhenoB
2     2  2000  2000 T     G     PhenoF,PhenoG,PhenoH,PhenoI
3     3  3000  3000 C     G     PhenoL,PhenoM,PhenoN

# The resulting table (tibble in this case) can be saved as a file
# using write.table() or write.xlsx()

```

`paste()` in combination with its `collapse` parameter is used because of the unequal numbers of columns in the files. `fread()` of the `data.table` package could be helpful for reading large files.
> 2 votes

# Answer
Convert your text files A and B to true, `sort-bed`-sorted BED and then use `bedmap` to map their associated IDs:

```
$ awk -vOFS="\t" '($1!~/#/){ $3+=1; print $1,$2,$3,$4; print $1,$2,$3,$5; print $1,$2,$3,$6; }' A.txt | sort-bed - > A.bed
$ awk -vOFS="\t" '($1!~/#/){ $3+=1; print $1,$2,$3,$4; print $1,$2,$3,$5; print $1,$2,$3,$6; }' B.txt | sort-bed - > B.bed

```

Then map:

```
$ bedops -u A.bed B.bed | bedmap --echo --echo-map-id-uniq --delim '\t' --multidelim '\t' - > answer.bed

```

If you need to put the header back in:

```
$ bedops -u A.bed B.bed | bedmap --echo --echo-map-id-uniq --delim '\t' --multidelim '\t' - | cat <(echo -e "#OUTPUT\n") <(awk '(NR==2)' A.bed) - > answer.txt

```

Some work might be needed with `awk` to put reference and alternate alleles back into their own columns, but this should get you about 90% of the way there.
> 1 votes

# Answer
**Edit:** I could not get bedtools to produce exactly the right output. See note at the end.

First install bedtools. Then, in Bash:

```
(tail -n +2 FileA.txt; tail -n +2 FileB.txt) | sort -k1,1V -k2,2n -k3,3n - > FileC.txt
head -n 1 FileA.txt > FileD.txt
bedtools merge -i FileC.txt -c 4,5,6 -o distinct,distinct,distinct -delim ", " >> FileD.txt

```

First line: `tail` takes lines 2+ of the file, `sort` sorts the files based on the chromosome, start, and stop coordinates, `-` is the output from the previous command.

Next line: `head` takes the first 1 line of FileA (header), saves to FileD.

Last line: `bedtools merge` merges the intervals that overlap. Option `-c` tells `merge` what columns to do operations on when two intervals are merged. Option `-o` tells merge what to do with those columns; in this case, it will make a comma-delimited list of distinct values. Option `-delim` defines the delimiter for listing distinct info in columns 4-6.

Output:

```
#Chr    Start   End Ref Alt Phenotypes
1   999     1001    A   T   Pheno A, Pheno B
2   1999    2001    T   G   Pheno F, Pheno G, Pheno H, Pheno I
3   2999    3001    C   G   Pheno L, Pheno M, Pheno N

```

**PROBLEM:** This code will not produce the desired output b/c bedtools expands the interval `1 1000 1000` to `1 999 1001`. BED files are used for multi-nucleotide intervals whereas the VCF standard is better suited for single nucleotide intervals and "ref/alt" variant information. Bedtools might handle VCFs with 1-nt intervals appropriately, but but evidently not BED files with 1-nt intervals.
> 1 votes

# Answer
If your files are small enough to fit in memory, here's one way, in Perl:

```
$ perl -lane 'if(/^#/){$h=$_; next}; $k=join("\t",@F[0..4]); $v=join(" ",@F[5..$#F]); $d{$k} ? ($d{$k} .= ", $v") : ($d{$k} = $v) ; }{ print $h;for(sort keys(%d)){print "$_\t$d{$_}"}  ' file*
#Chr Start End Ref Alt Phenotypes
1 1000 1000 A T Pheno A, Pheno B
2 2000 2000 T G Pheno F, Pheno G, Pheno H, Pheno I
3 3000 3000 C G Pheno L, Pheno M, Pheno N

```

Since that might be a little unclear, here's the same thing as a commented script and with informative variable names:

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

## Iterate over each line of each input file. This is done by the -n switch.
while (<>) {
  ## remove trailing newlines. This is done by -l
  chomp;
  ## If this line starts with a #
  if(/^#/){
    ## Save it as a header. We do this for each header line, so only
    ## the last will be saved. but since the files both have the same header,
    ## that's not a problem.
    $header = $_;
    ## Skip to the next line. This is to avoid processing the headers.
    next
  }
  ## Split the line on whitespace into the @fields array. This is
  ## done by the -a switch which creates the @F array used above.
  @fields = split(/\s+/);
  ## The $key variable holds the concatenated string of the first 5 fields.
  $key=join("\t",@fields[0..4]);
  ## The $values variable holds the rest of the line.
  $values=join(" ",@fields[5..$#fields]);
  ## If we've seen this $key before
  if($data{$key}){
    ## Append the current data to the value stored in the %data hash for this
    ## $key.
    $data{$key} .= ", $values";
  }
  ## If we haven't seen this key before
  else{
    ## Save the values for this key in the %data hash.
    $data{$key}=$values
  }
}
## print the header
print "$header\n";
## Now that all the data have been loaded, iterate over
## the sorted list of keys from the %data hash.
foreach my $key (sort keys(%data)){
  ##And print the output
  print "$key\t$data{$key}\n"
}

```

If you save that as `foo.pl`, you can run it like this:

```
$ perl foo.pl fileA fileB
#Chr Start End Ref Alt Phenotypes
1 1000 1000 A T Pheno A, Pheno B
2 2000 2000 T G Pheno F, Pheno G, Pheno H, Pheno I
3 3000 3000 C G Pheno L, Pheno M, Pheno N

```

---

Note that I am assuming the headers `#File A` and `#FileB` aren't actually part of your input and that you don't want `#OUTPUT` included in the output file. If this isn't true, let me know and I can modify this.
> 0 votes

---
Tags: sequence-annotation
---


# Question
Title: How to input data for DESeq2 from individual HTSeq count?
I am comparing the gene expression of 2 bacteria under 1 condition. I have now the count tables for 3 tech. replicates for each bacteria.

```
Bacteria1_1.count
Bacteria1_2.count
Bacteria1_3.count

```

...same for the other bacteria.

These files look like this:

```
gene1 10000
gene2 500
gene3 0
gene4 5000

```

I want to use DESeq2 for differential gene expression analysis. But I cannot figure out how to properly execute the DESeqDataSetFromHTSeqCount() command with this type of data.

Is there another intermediate step to add?

# Answer
In case still needed, here is a piece of code that I have used for the same purpose (includes a lot of sanity checks). This is highly likely an overkill as each count file would probably have the same rows (as a result of being generated from the same workflow) and thus can easily be merged column-wise.

```
count_file_dir <- "../input/count_files/"
count_file_names <- list.files(count_file_dir)
count_files <- paste0(count_file_dir, count_file_names)

# read and merge count files, while removing last 5 lines that correspond to summary stats
count_table_list <- mclapply(count_files, function(x) {fread(x)[1:(.N-5)]})
names(count_table_list) <- sub(".count", "", count_file_names)

# check the dimensions of count files
t(sapply(count_table_list, dim))

# check if the count files have the same gene ordering
# compare gene names of the first file with all
gene_names <- sapply(count_table_list, function(x) x$V1)
as.data.frame(apply(gene_names, 2, function(x) all(x == gene_names[,1])))

# add file name as a unique column name for gene counts
for(i in names(count_table_list)) {
  names(count_table_list[[i]])[2] <- i
}

# merge read counts and convert to a matrix
count_matrix <- Reduce(merge, count_table_list)
count_matrix_rownames <- count_matrix$V1
count_matrix$V1 <- NULL
count_matrix <- as.matrix(count_matrix)
rownames(count_matrix) <- count_matrix_rownames
colnames(count_matrix) <- sub(".count", "", colnames(count_matrix))

dim(count_matrix)

dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                              colData = ...,
                              design = ...)

```
> 2 votes

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


# Question
Title: Total read count of shotgun affected by few highly abundant sequences
I am building statistical models to analyse output from Illumina shotgun sequencing (HiSeq 4000) on stool samples (but RNA-seq data should behave similarily). The raw counts are a statistical sample of the total collection of DNA molecules used as input for the sequencing machine, and I think we can assume that the sequencing platform has limited capacity, such that detected reads are always a small fraction of the the original sequencing library.

Consider two different samples: one control sample and one treatment sample, which is identical apart from a single species (or gene) that has much higher absolute abundance relative what it in the control. As a result, the reads of this species will be much higher in the treatment sample, and also the 'real estate' of the remaining genes in that sample will be decreased. This is simply the nature of relative abundances.

However, size factor estimation that uses the "median ratio method" (http://dx.doi.org/10.1186/gb-2010-11-10-r106) for instance, are based on "experience with real data \[...\] \[that\] shows \[that\] a few highly and differentially expressed genes \[or abundant species\] may have strong influence on the total read count." (http://dx.doi.org/10.1186/gb-2010-11-10-r106). I am curious to see experimental evidence for this last statement. Or does someone have at least an explanation of how total read count can be affected by the absolute abundance of a species/gene?

# Answer
> 2 votes
That sentence in the Anders and Huber paper is a bit misleadingly worded (I don't think that was their intent). What they mean is that highly expressed and highly differentially expressed genes can out-compete the majority of the other genes, which thereby leads to them appearing to have lower expression if one uses total read counts for normalization (basically, what you said in the second paragraph).

An alternate interpretation of what they wrote (likely the one you understandable understood to have been their meaning) is that some highly expressed genes will inflate the total number of reads a sample actually gets. I have yet to see this in practice, but it certainly COULD happen. There are technical biases in sequencing such that some sequences (often due to length and GC content) are easier to sequence on common machines. A good example of this would be PhiX, which forms clusters much more easily on Illumina machines than real samples. I expect that this is relatively rare in practice, but using a more robust normalization would still help in protecting against it.

---
Tags: statistics, illumina, coverage, read-correction
---


title: Retrieve detailed gene descriptions

text:
# Question
Title: Retrieve detailed gene descriptions
Given a list of gene IDs, how do you retrieve the gene description, summary and other detailed information in R?

# Answer
## MyGene.info

Using the R package for the MyGene.info database:

```
source("https://bioconductor.org/biocLite.R")
biocLite("mygene")

library(mygene)
gene <- getGene("1017", fields="all")
gene$summary

```

Multiple genes:

```
gene <- getGenes(c(1347:1350), fields="summary")

```

---

## Entrez database

```
library(devtools)
install_github("ropensci/rentrez")

library(rentrez)
genes <- entrez_summary(db="gene", id=c(718,720))
genes$`720`$summary

```
> 7 votes

# Answer
Description of how to use biomaRt here. Filter by your list of gene IDs. Get the descriptions etc as attributes. Eg:

```
> library(biomaRt)
> ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
> IDs <- c("BRCA2","BRAF")
genedesc <- getBM(attributes=c('external_gene_name','description'), filters = 'external_gene_name', values = IDs, mart =ensembl)
> genedesc
  external_gene_name                                                                      description
1               BRAF B-Raf proto-oncogene, serine/threonine kinase [Source:HGNC Symbol;Acc:HGNC:1097]
2              BRCA2                  BRCA2, DNA repair associated [Source:HGNC Symbol;Acc:HGNC:1101]

```
> 8 votes

# Answer
The easiest and fastest method (which also doesn’t require internet access) is to use annotation packages, e.g. org.Mm.eg.db:

```
library(AnnotationDbi)
library(org.Mm.eg.db)

# 1:1 mapping symbol→entrez
entrez_ids <- mapIds(org.Mm.eg.db, gene_symbols, 'ENTREZ', 'SYMBOL')

# 1:n mapping symbol→…
print(columns(org.Mm.eg.db))  # check what you can query
df <- select(org.Mm.eg.db, gene_symbols, c(...), 'SYMBOL')

```
> 2 votes

---
Tags: r, gene, sequence-annotation, public-databases
---


title: Chromosome information for miRNA

text:
# Question
Title: Chromosome information for miRNA
Apologies if this question is previously asked. I have a `BAM` file alignment to `hg19` and miRNA gff3 file from `mirbase` for annotation of BAM file. I want to extract information in the following format (e.g.):

```
miRNA        chr    chr_start    chr_end    expression_count    sequence
hsa-mir-155  chr2       ...         ....         1020            ATGC....

```

I know that I can use htseq-count for BAM annotation and for expression counts but how can I get their chromosomal location. The reason for doing this specific way is that I want to show them in a CIRCOS plot.

Just counting is not the goal, their chromosomal location also needs to be extracted. So just subread/htseq-count is not sufficient for extracting genomic location information.

Previously, I had the same idea about processing the sam records. If I want to annotate chromosomal information with miRNA using .gff file, I have to check for each line in .sam records whether their chr\_start and chr\_end lie in the bandwidth of start and end of any miRNA in .gff file. This process can be very heavy and slow because .sam records are of more than 5 GB. That's why I want to know if there is any extremely fast and generalized way to do it.

I want to avoid processing sam records in python because my experience with extremely heavy data like fastq or sam in python is not good. I would rather prefer awk/bash or any low-level language.

Thanks.

# Answer
you can do `samtools view pathtobamfile`, this will convert the file to `.sam` format and send it to `stdout`, the columns can then be parsed out which contain the information you are after.
> 1 votes

# Answer
The gff should have the coordinates for each feature. I'd use featureCounts from subread to count up how many reads hit each feature.
> -1 votes

---
Tags: sam, circos, chromosome-capture
---


# Question
Title: Compare multiple alignment results' aligned bases
I have aligned a nanopore data set to a reference genome with graphmap, minimap2 and BLASR. The alignment results are stored in BAM files.

I would like to do some concordance assessment, looking at the number of base pairs that are mapped to the same position by different aligner.

Are there any tools or ways that allow me to do this effectively?

# Answer
> 1 votes
I don't know of a single tool that can read multiple sets of BAM files and compare them. But you can use the Sam\_Reader python3 module to create a csv file that summarizes each position in the genome for each of your 3 results and then compare them manually.

It should be installable with:

```
pip3 install sam_reader

```

And then generate the CSV with

```
from sam_reader import Sam_Reader
my_files = Sam_Reader('my_bam_file_or_folder/')
my_files.per_base_stats(write_file='new_file.csv')

```

\*but FYI this module requires pysam, which does not work on windows

# Answer
> 1 votes
I'm not aware of any tools that do this out-of-the-box. Assuming some programming experience, it would be fairly straightforward to implement this in Python using pysam (or htslib bindings for your favorite scripting language). In particular, the pileup is a very useful data structure.

Ultimately, you'll need to store an index of reference sequence positions by read name and position.

```
read 1, offset 0 --> chrom 3, offset 248879
read 1, offset 1 --> chrom 3, offset 248880
read 1, offset 2 --> chrom 3, offset 248882
...
...
...

```

Once you've done this for a pair of BAM files, computing the concordance statistics involves trivial set operations.

There are some important practical considerations though. If you implement the approach completely in memory, it will only be tractable for small-ish BAM files. It wouldn't take very many reads (relatively speaking) to overwhelm available RAM on a laptop or desktop. For large BAM files, you might even overwhelm the available RAM on a big-memory server. Alternatively, you could store the indexes in an on-disk database such as sqlite3. This would keep RAM requirements low, but would slow things down substantially.

I think the most reasonable approach would be to process the BAM files in batches: grab all the reads in a 100kb region (or whatever size makes sense), compute the index, store the congruence statistics, move on to the next 100kb region and repeat, and then aggregate statistics when all regions have been processed.

# Answer
> 0 votes
You can use the bam\_alignments\_compare.py script from the wub Python package.

---
Tags: sequence-alignment, sam, nanopore
---


# Question
Title: Concatenate multiple fasta files into one file per unique entry
I have eight fasta files, one per gene, with 106 fasta entries per file, like this:

gene1.fasta:

```
>sample1_434564
GATCTCGTAGATGTGA
>sample2_2346764
GATCTCGTAGATGTGA

```

gene2.fasta:

```
>sample1_532345
ATGCTGTGAAGTA
>sample2_6040103
ATGCTGTGAAGTA

```

What I want to do is create one fasta file per sample so that I get all gene sequences for that sample in one file. As you can see, the sample ID is the same before the `_` in the fasta headers, but the rest is unique for each fasta file.

Is there any way to do this with `awk` or other command-line tools?

Desired output:

Sample1.fasta

```
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

```

Sample2.fasta

```
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

```

# Answer
You can do this simply enough in `awk`:

```
$ awk -F'_' '{
               split(FILENAME,n,/\./);
               if(/^>/){
                  sample=$1;
                  sub(/>/,"",sample);
                  $0=">"n[1]
               }
               print > sample".fasta"
             }' gene1.fasta gene2.fasta ... geneN.fasta

```

If I run this on the two files you gave, I get:

```
$ cat sample1.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

$ cat sample2.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

```

### Explanation

* `-F'_'` : set the field separator to `_` as an easy hack that lets us use `$1` for the sample name (`sample_NNN`, so `sample` will be `$1` and `NNN` will be `$2`).
* `split(FILENAME,n,/\./);` : `FILENAME` always holds the name of the file currently being processed. Here, we are splitting it on `.` to get the gene's name. **This assumes all of your files are named `$geneName.fasta`.** The `split` creates the array `n` whose 1st element, `n[1]` is the gene name.
* `if(/^>/){` : if this is a fasta header line.
* `sample=$1` : set the variable `sample` to the sample's name. See the first note about `-F'_'`.
* `sub(/>/,"",sample);` : remove the `>` from the `sample` variable. We don't want that in the file name.
* `$0=">"n[1]` : set the current line to be `>` and the gene name.
* `print > sample".fasta"` : print each (modified) line into the appropriate file name.
> 6 votes

---
Tags: fasta, text-processing
---


title: Merging ssu-align alignments

text:
# Question
Title: Merging ssu-align alignments
ssu-align is an HMM-based MSA tool for ribosome's small subunit (SSU) sequences. It's bundled with three SSU models:

* archaeal,
* bacterial
* eukaryotic

**Approach** Given a set of archaeal and bacterial sequences the tool will first identify the best model for each sequence and produce two separate alignments for archaea and bacteria, and the tool provides no merging functionality.

**Rationale** The underlying issue with SSU is the secondary structure, which results in: 1) indels; 2) pairwise point mutations (to maintain the secondary structure)

**Problem** While you can force it to use a single model for all sequences, aligning archaeal sequences to the bacterial model will produce poor results.

**Solution?** Nevertheless, the authors of GreenGenes claim that they've used ssu-align to produce two MSAs (one for archaea and bacteria) and build a joint phylogenetic tree using FastTree, but they say nothing about the way they've managed to merge the alignments. Am I missing some functionality from SSU-align?

**Note** Merging two giant MSAs using conventional tools, such as MAFFT, won't get me far.

P.S.

I've already tried contacting the authors of the GreenGenes paper, but the corresponding author hasn't responded.

# Answer
I think this is more to do with tree building than alignment construction.

I believe the underlying issue is the 1) distribution of indels, this results in a blockwise deletion in the resulting phylogeny and 2) nucleotide saturation.

**Indels** Thus, for indels

```
Archaeal
XXXX----YYYYY
bacterial
----ZZZZYYYYY

```

The resulting phylogeny will be based on YYYYY and XXXX is ignored.

The alternative explanation is that the genetic distances involved have resulted in saturation of the nucleotide sequences. SSU has to be a nucleotide phylogeny and the level of reversion mutations over a vast amount of evolutionary time can become indistinguishable from random.

With regards the first point, you can replace the indels "-" with "N". This is a cheat to trick a tree building program into not discarding the indel. Model based tree building programs have no ability to create a genetic distance between a base pair and an indel, so this data is essentially discarded from the tree. Worse it will remove the site from phylogenetic inference - which can be a massive loss of data. Two notes one this,

* Some max. likelihood, Bayesian tree programs can avoid blockwise deletion, I don't know programmatically how they do this;
* In PAML the distribution of indels was modelled, it doesn't appear to have been widely adopted however, but might be of value.

**Alternative** The alternative method is really for protein genes and called "RY" coding, it might work on SSU because of the compensatory mutations involved in maintaining a secondary structure of the SSU.

**Word of warning**, the saturation may be so great, very little nucleotide point mutation is preserved. Essentially we are talking about preservation of phylogenetic signal across the origins of life ... and that is a very, very long time ago.
> 0 votes

---
Tags: sequence-alignment, fasta, phylogenetics, phylogeny, rna-alignment
---


# Question
Title: Can I change the name of file features.tsv to genes.tsv
it is said that the features.tsv from Cell Ranger v3 is analogous to the genes.tsv from Cell Ranger v2. So can I change the file name to genes.tsv for Seurat to read it? I know that Seurat has the update version 3 that can read the features.tsv, and I would like to use the old version. So if I change the file name, did it cause any problem?

# Answer
> 1 votes
You can follow the 10X instructions to load your matrix into R: the following is a direct quote from their webpage:

> ```
> library(Matrix)
> matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
> barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
> features.path <- paste0(matrix_dir, "features.tsv.gz")
> matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")
> mat <- readMM(file = matrix.path)
> feature.names = read.delim(features.path,
>                            header = FALSE,
>                            stringsAsFactors = FALSE)
> barcode.names = read.delim(barcode.path,
>                            header = FALSE,
>                            stringsAsFactors = FALSE)
> colnames(mat) = barcode.names$V1
> rownames(mat) = feature.names$V1
>
> ```

You can then create a seurat object using mat

```
object <- CreateSeuratObject(mat, project="")

```

---
Tags: single-cell
---


# Question
Title: How to handle control samples in CLIP-Seq
I have a CLIP-Seq dataset I'm processing, which includes control samples and no inputs. This is the second CLIP analysis I've performed to help out users of our genome core facility and the first one only had inputs. I've done a number of ChIP analyses in the past, which I loosely consider to be a similar type of analysis, and they all used inputs, so I've never had to look into the differences in using inputs versus controls when doing peak calling before. Now that I'm trying to figure out what I need to do, I'm a little lost.

Here's what I've got:

1. Samples of cross-linked RNA/protein complex that went through a poly-A and then a FLAG pull-down (protein G beads with an antibody). The complex was eluted, treated with proteinase K, and the RNA was sequenced.
2. Controls went through the same processing as the samples except the control beads had no antibody

There are replicates in both cases. (As a side-note, the RNA was treated with RNase III while bound to the beads and the RNA's trimmed off the complexes were also sequenced. I'm not sure what to do with those samples yet.)

I was going to use one of the CLIP peak callers available on Galaxy, but whichever one I use, all the protocols I've found and the documentation I've looked at use inputs for the peak calling method and not controls. So what do I need to do differently (if anything) to use the controls as the background in the peak calling? Or do I need to suggest they also sequence the inputs?

And just so I know (dumb question here)... Is "input" (as the name suggests) the sample as it was before the FLAG pulldown? Why would one use one or the other (input versus control).

# Answer
I've not handled CLIP-seq before, but we have done iCLIP and in this we had controls, but no input. @DevonRyan is correct, these control samples are similar to IgG pulldowns in ChIP seq.

The way we handled this was to call our peaks on isolated samples. The original peak calling alog from the first CLIP seq paper doesn't require input. Also I believe that clipper doesn't need input (although it will take it). We then removed any peaks from the sample that overlapped a peak in the control.

But most of our analyses we didn't do on peaks - we did lots of things like looking at metagenes, or classes of genes, or particular types of exon etc, and when we did this, because we were averaging over many sites, we could generate the pattern for the control as well as for the input, and often normalise the one to the other.
> 0 votes

---
Tags: rna, peak-calling
---


# Question
Title: Are there any resources for calculating / representing evolutionary relationships within gene clusters?
I am working with a gene cluster found in several species. The genes involved are pretty much the same, with the odd insertion, but they can be often found in different orientations. So, I can't just do a multiple alignment of the whole cluster. I've done phylogenies of individual genes and I've looked at cluster architecture (indels, rearrangements) and things make sense (eg. limited HGT), but..

...Is there a simpler way to do this? A way to calculate / represent a phylogeny of the whole cluster? Like...some sort of web, maybe? Most of the software that does take rearrangements into consideration is meant for whole genomes, and is unreliable when it comes to fragments about 20-30kb long (like mine).

Thank you for your help!

# Answer
Firstly you wanted a webserver to build a tree, ok, I would use PhyML server. The good thing is that with only 22 taxa the calculations are easily doable.

What I'm assuming is you can map your rearrangements to produce a single homologous gene alignment for all taxa. So you appear to infer you can break your 20-30kb contig into individual aligments of homologous genes which has incorporated the rearrangement and enables easy production of an alignment file. This would leave you with around 20-30 gene alignment files. I hope that is correct, because that is how you get around the rearrangement problem and enables the data to be reassembled for the final tree (below).

Doing the tree is simple, because you can exploit the assumptions used in tree building regarding independence between mutations at each alignment position. The theoretical proof, the data is sound, is complicated. In summary,

* build a consensus phylogeny (one tree representative of lots of trees) using two different approaches
* combine the data into a single alignment, if the above looks good, or excluding the alignments that are affected recombination.

The theoretical issues with what you want to do are complicated because you need to distinguish bifurcating (clonal) homologues from recombination that could have occured in the reassortment events. The following is a purists approach, you might not want to do this stuff, but it gives you some idea.

Just before I do that ... I just want to check the alignments assumption is correct? and which bacteria?
> 1 votes

---
Tags: fasta, phylogenetics, gene, phylogeny, pathway
---


title: Comparing RNA-Seq and Ribo-Seq

text:
# Question
Title: Comparing RNA-Seq and Ribo-Seq
I'm attempting to compare total RNA-seq with Ribo-Seq, to determine if changes in Ribo-Seq are due to changes in transcriptional expression (akin to the analysis performed by Anota2Seq). I am, however, struggling to understand how to accomplish this in DESeq2.

My design is as follows:

```
cts <- read.delim("Interaction.txt", row.names=1)
coldata <- read.delim("Sample_Sheet.txt", row.names=1)
cts <- round(cts)
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)
dds <- DESeq(dds)
resultsNames(dds)

> coldata
                Genotype Type
C_24h_Ribo_1     Control Ribo
C_24h_Ribo_2     Control Ribo
C_24h_Ribo_3     Control Ribo
CHIK_24h_Ribo_1    Viral Ribo
CHIK_24h_Ribo_2    Viral Ribo
CHIK_24h_Ribo_3    Viral Ribo
C_24h_RNA_1      Control  RNA
C_24h_RNA_2      Control  RNA
C_24h_RNA_3      Control  RNA
CHIK_24h_RNA_1     Viral  RNA
CHIK_24h_RNA_2     Viral  RNA
CHIK_24h_RNA_3     Viral  RNA

> resultsNames(dds)
[1] "Intercept"                 "Genotype_Viral_vs_Control" "Type_RNA_vs_Ribo"          "GenotypeViral.TypeRNA"

```

I do not understand which contrast is the one I need, or how to access the results i'm looking for (or interpret them). Any help?

# Answer
Adding to what benn pointed out, at this line:

```
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)

```

You specified that the counts, can be a explained by Genotype (being Control or Viral) or Type (RNAseq or Ribo) or an interaction effect

What you need is this coefficient, Genotype:Type. To use it, go

```
res <- results(dds,name="GenotypeViral.TypeRNA")

```

Based on your setup, you are always comparing RNA to Ribo (see "Type\_RNA\_vs\_Ribo" in resultsNames). If changes to the abundance of a mRNA is similar in RNAseq or Ribo, then the coefficient of the interaction term will be zero and it will be zero in the log2foldchange column in results(res). A positive fold change indicates the log2foldchange between RNA and Ribo is higher in Viral.

If you would like to compare ribo versus and rna-seq, then set RNA as the reference by:

```
coldata$Type = factor(coldata$Type,levels=c("RNA","Ribo"))
# run deseq again
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)
dds <- DESeq(dds)

```

With this you will compare by default Ribo versus RNA, and the interaction term will tell you how much more enriched or depleted the ribo fraction is, in the Viral samples.

Would recommend checking the vignette out here:

http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts
> 3 votes

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


title: Mapping probe ID to gene name

text:
# Question
Title: Mapping probe ID to gene name
I have set of modules from various condition ,with first column as Probe\_ID name and I have an annotation file with Probe\_ID and Gene\_name.

I have to map the Probe\_ID to gene name in each file and write it back.

So what im doing is as such

```
filelist = list.files(pattern = ".*.txt")

```

I get list of files

Then those list of files into list of dataframe

```
datalist = lapply(filelist, function(x)read.table(x, header=T))

```

Now i read into the annotation file

```
ANNOTATION_FILE <- read_delim("/run/media/punit/data1/Dev_bio_info/ANNOTATION_FILE.txt",
    "\t", escape_double = FALSE, trim_ws = TRUE)

```

The normal way i do is

I perform inner\_join(ANNOTATION\_FILE,my\_file) for single file but as i have like 22 files in the module .I want to do it in loop and im not sure how to proceed further.

My data frame example

This is how my module files looks like .

```
dput(head(module_ivory[c(1:10)]))
structure(list(Gene = c("1556672_a_at", "1557803_at", "1558080_s_at",
"1566465_at", "213089_at", "213229_at"), GSM1304905 = c(5.977158746759,
6.83852532860115, 5.93931007942615, 6.11926533584282, 9.92182534751875,
10.2576150019626), GSM1304906 = c(7.51622774902859, 7.71268478229187,
7.32360342148457, 6.8928342169826, 9.12807784837542, 10.6767090567628
), GSM1304907 = c(6.58069278884036, 7.86890363253932, 6.9901005950279,
6.08614452085924, 10.1220627669628, 10.9203125158285), GSM1304908 = c(6.60119000844375,
6.71735654004629, 7.43132441110608, 5.48084261893207, 9.81689665710992,
10.7649566455486), GSM1304909 = c(5.95049879217177, 6.132975671606,
7.0902951305213, 5.59408027736444, 9.08474568290096, 9.84784203563294
), GSM1304910 = c(5.83616005353191, 5.94947517647125, 6.61469427098759,
4.31285839654272, 9.44446310110043, 9.99021299504538), GSM1304911 = c(5.92373056182016,
6.14502105445603, 6.2333758614231, 5.46455403294761, 9.55447963883665,
10.1810888944642), GSM1304912 = c(6.43107367158601, 6.36340779490876,
7.18374247647159, 5.71508717609646, 9.61355221684561, 10.3266656836978
), GSM1304913 = c(6.4010171640209, 7.31600259413326, 6.16969520557938,
5.40234538288784, 8.89481478840621, 10.290064262843)), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -6L))

```

My annotation file

```
dput(head(ANNOTATION_FILE))
structure(list(Gene = c("1007_s_at", "1053_at", "117_at", "121_at",
"1255_g_at", "1294_at"), Gene_Symbol = c("DDR1 /// MIR4640",
"RFC2", "HSPA6", "PAX8", "GUCA1A", "MIR5193 /// UBA7")), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -6L))

```

I would like to do inner join on each file and write them back into the file .

Any help or suggestion would be really appreciated

# Answer
> 1 votes
To end up with a single table made by merging all your tables, you can give a try with `reduce()` of `purrr`: The first argument `.x` would be the list of tables to merge (in your case the annotation file plus all the other files) and the `.y` argument would be your function, i.e. `inner_join()`. `reduce()` would merge all the tables within the given list, left to right or right to left. `...` would enable you to specify column(s) to use as "anchors" while merging.

If you would prefer to keep your tables separately then `lapply()` is the answer: `lapply(list_of_files, function(x) inner_join(annotation_file, x))`. This would give a list of merged tables.

---
Tags: r
---


# Question
Title: Program to make a haplotype network for a specific gene
Currently my lab manually uses pop art to make the network after using PHASE to id an individual's alleles. I would like to know how to generate the graph structure (i.e. the links between vertices), as well as how to draw the graph.

Input data are formatted as FASTA files with the sequence in nucleotide formatting.

Example:   I just feel like there's got to be a better, more time effective way to do this.

# Answer
> 1 votes
This is a haplotype map, where each node is proportional to the frequency of a given allele. The formal tool to do this is Network 5 found here. The problem with manaul haplotype maps is the area of the node should be proportional to the frequency, if this is not accurate it could be perceived to be misleading.

The approach for looking at each parsimonious state along a branch is a phylogenetics method which is performed via MacClade 4.0 ... requires an old Mac though.

---
Tags: networks, phylogeny, phylogenetics, haplotypes
---


# Question
Title: How to find common sequences among 6 multi-fasta files
I have **6 multi-fasta files**, every of them contains ca 1500 sequences like that:

```
>Haladaptatus sp. R4
MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
ALEEVADR
>Methanosphaera sp. WGK6
MAKVNVYSLKGDITEEIELPEIFEEEYRPDVIKRAVISTQTARIQPWGANPMAGKRTTAE
SFGSGRGAAMVPRVKGSSKAAFVPQAIGGRKAHPPRVNTIYHEKINKKERILAIRSAIAA
TANKEIVEQRGHAVANLEQVPFVVDDELETIKTTKETREIFKDLGIMDDILRAKKGRKIK
SGKGKLRGRKYRTPKGPLVVVGNDRGISLGARNHAGVEVVEVNNINAELLAPGTHAGRLT
IYTKSAVEKLADLFQQNRS

```

and so on. I need to create a new result file containing sequences that are common in these multi-fasta files. Example: If sequence

```
>Haladaptatus sp. R4
    MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
    PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
    AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
    GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
    ALEEVADR

```

occures in 1 **and** 2 **and** 3 **and** 4 **and** 5 **and** 6 file it shoud be written to the result file. May anoyone help?

Edited: this is my Python script which should work on 2 files (not 6 for now to make it simplier):

```
for z in f1:
if z[0] == '>':
    znaleziony = False
    gatunek1 = z[1:len(z)]
    #print(gatunek1)
    for x in g1:
        #print(x)
        if x[0] == '>':

            gatunek2 = x[1:len(x)]
            if gatunek1 == gatunek2:
                wynik.write('>')
                wynik.write(gatunek1)
                znaleziony = True
                continue

elif znaleziony == True:
    wynik.write(z)

```

wynik.close()

Unfortunately it repeats some identifiers of sequence and does compare only identifiers (not whole sequence).

# Answer
> 3 votes
If both the ID lines and the sequences are duplicated, you can do this using the FastaToTbl and TblToFasta scripts I have posted before and just normal Linux utilities:

```
FastaToTbl file1.fa file2.fa ... file6.fa |
    sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //' | TblToFasta > dupes.fa

```

To understand this, let's look at the output of each step. The `FastaToTbl` script just converts fasta to a tab-delimited format where the ID is the first field and the sequence is the second, all on one line:

```
$ FastaToTbl file*fa | head -n1
Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

```

Then, we pass this through `sort` and `uniq -dc` which prints out the number of times a line was repeated (only for duplicated lines):

```
$ FastaToTbl file*fa | sort | uniq -dc
      4 some other sequence ALASOOPAKASKAPPPQQCDDREETAYSSLDAS
      6 Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

```

Since we only want those that were present in all 6 files, we can use `awk` telling it to only print lines whose first field was 6:

```
$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6'
      6 Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

```

But then we don't want the 6 as part of the final output, so we need to remove it:

```
$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //'
Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

```

Finally, convert to fasta again:

```
$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //' | TblToFasta
>Haladaptatus sp. R4
MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
ALEEVADR

```

---

Alternatively, if your files are small enough to be loaded into memory, you could do this instead:

```
FastaToTbl file*fa |
    awk '{a[$0]++}END{for(s in a){if(a[s]==6){print s}}}' |
        TblToFasta > dupes.fa

```

Or you could just do the whole thing in a single awk operation:

```
awk '{
        if(/^>/){
            a[seq]++;
            seq=$0;
        }
        else{
            seq=seq"\n"$0
        }
      }
      END{
        a[seq]++;
        for(s in a){
            if(a[s]==6){
                print s
            }
        }
       }' file*.fa

```

# Answer
> 0 votes
For pandas you simply load the fasta names and fasta sequences into separate arrays. The "fasta names" array becomes the index in a pandas Series, whilst the sequence are the values. The dataframe will correctly re-associated the names with each sequence. You load all the files into a single dataframe and use the duplicated() function then write out.

The easiest way to do this is concatenate all your files into one via shell ...

```
cat file*.fasta > concat_file.fa

```

Python script (biopython + pandas) ...

```
from Bio import SeqIO
import pandas as pd

path = "/Users/username/location/concat_file.fa"
outpath = "/Users/username/Desktop/duplicates.fasta"
fasta_name = []
fasta_seq = []

for record in SeqIO.parse(path, "fasta"):
    fasta_name.append('>' + record.id)
    fasta_seq.append(record.seq)
df = pd.Series(fasta_seq, index = fasta_name)
df_dups = df.index.duplicated()
mydups = df[df_dups]
mydups.to_csv(outpath, sep="\n", index=True, header=False)

```

...

I checked it on a ~4000 sequence fasta file of unique sequence ids (dengue), duplicated the last sequence id + sequence and the only sequence output was the very last sequence.

---
Tags: python, fasta, ngs, phylogenetics, phylogeny
---


# Question
Title: After running a DRIMSeq pipeline, how do I know which genes are upregulated in the different conditions?
After running a DRIMSeq pipeline and obtaining the genes that are differently used between the null and full models, how do I select the genes that are differently used in the different conditions?

According to the vignette, I would have to use the `plotProportions` function and say in the `group_variable` the variable of my interest. What if I want a list of the genes that are differently used specifically in the variables of my interest? Checking every gene manually sounds exhausting. Help would be much appreciated.

# Answer
> 0 votes
After running the test for differential transcript usage you get a list with gene- and transcript-level p-values that you can use to filter which differences are significant. The gene-level p-value tells you if there is differential usage of transcripts within a gene and the transcript-level p-value tells you wich transcripts actually are changing. Therefor you will want to select significant genes for further investigation.

The plots are mainly used to visualize top hits as detailed in the vignette.

```
res <- results(d)
# The data is ordered by p-value
res <- res[order(res$pvalue, decreasing = FALSE), ]
# The top hit is selected
top_gene_id <- res$gene_id[1]
plotProportions(d, gene_id = top_gene_id, group_variable = "group")

```

However you should consider running the two-stage test that is described on page 16 of the vignette you are using to reduce the number of false positive hits. From this you will get another list of adjusted p-values that can be used as above.

---
Tags: rna-seq, bioconductor
---


# Question
Title: snakemake get environment variable for rule name and log path
When a rule fails in snakemake workflow, it returns a descriptive message in log file like following:

```
Error in rule step1:
    jobid: 99
    output: results/foo/step1/sample5/outdir/output1.tsv, results/foo/step1/sample5/outdir/output2.tsv
    log: logs/step1/sample5_out.log (check log file(s) for error message)
    shell:
    ...
        (exited with non-zero exit code)

```

Since I run workflow over many samples, I have defined a *postrun.sh* bash script in *jobscript.sh* directive that runs after snakemake `{exec_job}` section and notifies me via userbot in slack channel when a step failed for one or more samples.

**jobscript.sh**

```
        ########################### START SNAKEMAKE CMD ############################
    {exec_job}
        ############################ END SNAKEMAKE CMD #############################
    exitstat=$?
    export exitstat

    echo "END at $(date)"
. "${{HOME}}"/.smk_confs/postrun.sh

```

**postrun.sh**

```
## execute after running each of snakemake job
echo 'Running postrun.sh'
ERRMSG="MYJOB ID: $PBS_JOBID failed at $(pwd) on $(hostname) for ${USER} with exit status: ${exitstat}."

## ping to slack
ping_slack -m "$ERRMSG"
sleep 5

```

I am looking for a snakemake exported environment variable that I can use in *postrun.sh* to know which rule ID failed for one or more samples. For now, I can think of using job scheduler variable matching job name, e.g, PBS\_JOBNAME, but better if I have a variable from snakemake API for rule ID or preferably, path to log file.

Thanks

# Answer
> 1 votes
You could add an exit trap to the beginning of the shell scripts inside your Snakefile, which is executed after each script is done and sends a message when there was an error:

```
 shell:
    """
    trap 'if [ $? -neq 0 ]; then ping_slack -m "Error in rule xxx with input {input}"; exit 1' EXIT
    ...

```

---
Tags: snakemake
---


# Question
Title: How does one construct a cladogram of intraspecies relationships?
I have SNP data from several cultivars of rice which I have used to produce alignments, but I don't think that the usual models and algorithms used for generating phylogenetic trees are appropriate, because these cultivars are not the result of speciation events and have been interbred in their histories. How can I best calculate and visualize their degree of relatedness?

# Answer
> 3 votes
Even with inbreeding and other genetic phenomena that might mask actual evolution of these cultivars, any phylogenetic methodology would be capable of determining relationships accurately.

Try creating a Neighbour Joining tree with MEGA, which is one of the simplest methods available. This should give you enough to check the relationships of the cultivars.

# Answer
> 1 votes
If the accessions have been interbred (conflicting with a tree-like phylogeny), you could try Splitstree, which tries to estimate a tree-like phylogeny with added hybridisation events.

# Answer
> 0 votes
Wait, you are right, but trees don't work here.

This is the area of population genetics. You will generate allele frequencies for each locus. The allele frequenceies will be used to perform Hardy-Weinberg equilibrium, F statistics, Fst and Fis

H-W start by using Structure API to work out how many populations you have. This is vital information for the next two analyses:

* Fis is particularly useful here because it is the inbreeding co-efficient, the nearer 1 it is the more inbreeding the data has been through.
* The Fst distances can be used to generate a tree, but this is not a model of point mutations, rather a model of migration between populations. You can input the distance matrix in MEGA and make a tree from it, Fstat API (I think) will do it for you.

Pop gen is difficult to get your head around at first, but you'll get use to it.

---
Tags: fasta, phylogenetics, snp, phylogeny
---


title: PDB codes for misfolded proteins?

text:
# Question
Title: PDB codes for misfolded proteins?
I'm searching for PDB code examples for "healthy proteins" and misfolded ones, so I can compare the 3D structure. I'd prefer pairs of a "protein misfolding disease", but could be any. Can I find such a list or how should I search for the particular codes? I'm particularly interested in ones that also cointain SS-bonds.

# Answer
> 1 votes
PDB `5E83` sickle cell anaemia.

All details for comparison to wildtype rbc are found in,

Oksenberg, D., Dufu, K., Patel, M.P., Chuang, C., Li, Z., Xu, Q., Silva-Garcia, A., Zhou, C., Hutchaleelaha, A., Patskovska, L., Patskovsky, Y., Almo, S.C., Sinha, U., Metcalf, B.W., Archer, D.R. GBT440 increases haemoglobin oxygen affinity, reduces sickling and prolongs RBC half-life in a murine model of sickle cell disease. (2016) Br.J.Haematol. 175: 141-153

# Answer
> 2 votes
Searching "misfolded" on the RCSB PDB turns up 64 hits. Most of these won't be relevant, but there are some interesting entries on tumor associated EGFR and misolded RNA.

In order to experimentally obtain a structure of a protein it generally has to be in a stable conformation, i.e. at low energy. Being soluble also helps. This usually corresponds to the "healthy" form, as misfolding would lead to aggregation.

That said, the structures of certain aggregates have been obtained. For example, see this recent article.

There is also the concept of mutating a protein sequence to fold to a different structure, see for instance here, though I guess you are looking for cases with the exact same sequence.

# Answer
> 2 votes
There are very few real structures of soluble misfolded proteins, for the reasons described by jgreener. An example with a disulfide bond is beta-2-microglobulin, eg PDB entries 2XKS vs 2XKU.

Atomic-resolution structures of amyloid fibrils may be what you're looking for, such as 6GK3 for beta-2-microglobulin amyloid (which has its disulfide intact). A search for "amyloid fibril" at rcsb.org will find many structures. But many amyloid precursor proteins are intrinsically disordered, or are fragments of native proteins.

There's no list that I know of, but here's a recent review: https://www.ncbi.nlm.nih.gov/pubmed/30237470

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


title: samtools mpileup skipping read

text:
# Question
Title: samtools mpileup skipping read
I run the command:

```
samtools mpileup -O -s -q20-B -Q20 -f hg19.fa -r chr1:569929-569931 myFile.bam

```

and get:

```
chr1    569929  G   7   ...,,., EEEEEEE NTTVTTU 53,48,42,60,30,29,27
chr1    569930  C   6   ...,.,  EEAEEE  NTTTTU  54,49,43,31,30,28
chr1    569931  G   6   ...,.,  EEEEEE  NTTTTU  55,50,44,32,31,29

```

I am interested in **position 569930**.

I inspect the bam with:

```
samtools view myFile.bam chr1:569929-569931

```

and I see:

```
NS500355:NS500355:HHYN5BGXB:2:21311:2303:9255   99  chr1    569869  0   76M =   569900  107 CTCCATAACGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE<EEEEEEEAEEEEEEE    NM:i:2  MD:Z:5C3C66 MC:Z:76M    AS:i:66 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9321,76M,1;
NS500355:NS500355:HHYN5BGXB:4:22603:25857:2799  99  chr1    569877  45  76M =   569896  95  CGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCAC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE    NM:i:1  MD:Z:1C74   MC:Z:76M    AS:i:74 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9329,76M,1;
NS500355:NS500355:HHYN5BGXB:1:12105:11868:12011 99  chr1    569882  51  76M =   569900  94  CTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACC    /AAAAEEEEEAEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9334,76M,1;
NS500355:NS500355:HHYN5BGXB:2:13301:2681:6114   99  chr1    569888  9   76M =   569952  140 CTAGGCATACTAACCAACACACTAACAATATACCAATGATGGAGCGATGTAACACGAGAAAGCACATACCAACGCC    AAAA/E/EE//AAEEEEE/EAAEAA//EEAEEE/E/EEEE/////A/EE/EE/EE/E//EE/EEEE/E/AEA####    NM:i:4  MD:Z:6C19C15C29G3   MC:Z:76M    AS:i:57 XS:i:52 RG:Z:../H2O XA:Z:chrM,+9340,76M,5;
NS500355:NS500355:HHYN5BGXB:2:23110:10125:9334  99  chr1    569888  51  76M =   569904  92  CTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCC    AAA/AEE6EEEE6E///E/A/AAE/E/AEEEEEEE/EE/EEEAEEEEAEE//EEEEEEAEEEAEEAE/AEE/AEE<    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9340,76M,1;
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE<A/EE//AE/AA///AAE///EEEEE//6///E//EE/6EE/EEEEEEEA/A/A/E/E//AEAEEEEE6A6AAA    NM:i:0  MD:Z:58 MC:Z:59M17S AS:i:58 XS:i:53 RG:Z:../H2O XA:Z:chrM,-9340,18S58M,1;
NS500355:NS500355:HHYN5BGXB:1:21304:26480:15673 83  chr1    569890  11  76M =   569796  -170    AGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCAC    AEEEEEEEAAEEEEEEEEEEEE/EEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,-9342,76M,1;
NS500355:NS500355:HHYN5BGXB:4:11602:21237:3028  83  chr1    569900  51  76M =   569900  -76 ACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCAC    EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,-9352,76M,1;
NS500355:NS500355:HHYN5BGXB:3:13504:4369:8812   99  chr1    569901  51  76M =   569901  76  CCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9353,76M,1;
NS500355:NS500355:HHYN5BGXB:2:11203:19774:18704 99  chr1    569904  11  76M =   569952  124 ACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGT    6AAAAEEEEEEEEEEEEEAEEEEEEEEEEEEEEEE6EEEEAEEEAEEEEEEEEEEEEEEEEEEAE/E6EEAEEEEA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9356,76M,1;
NS500355:NS500355:HHYN5BGXB:2:11309:13954:17567 99  chr1    569905  11  76M =   569972  143 CACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9357,76M,1;
NS500355:NS500355:HHYN5BGXB:1:21202:13252:14049 99  chr1    569905  11  76M =   569927  98  CACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTC    AAAAAEE6EEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEEEE/EEEE<EEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9357,76M,1;
NS500355:NS500355:HHYN5BGXB:1:22311:21964:12072 83  chr1    569908  52  5S71M   =   569908  -71 GATCTACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTG    <AEEEA6EEE6EEEEEEEEEEEEEE/EEEEEEEEEEEEEEEEE/EEEEEEEEEEEEEEEEEEAAEEEAEEEAA6AA    NM:i:0  MD:Z:71 MC:Z:71M5S  AS:i:71 XS:i:66 RG:Z:../H2O XA:Z:chrM,-9360,5S71M,1;
NS500355:NS500355:HHYN5BGXB:4:12602:7858:19418  99  chr1    569909  15  76M =   569917  84  CTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTCCAAA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEE/AAEEEEEEEEEEEEEEE    NM:i:1  MD:Z:74G1   MC:Z:76M    AS:i:74 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9361,76M,1;
NS500355:NS500355:HHYN5BGXB:2:12205:19456:7598  83  chr1    569927  0   76M =   569907  -96 TGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTCCAAAAAGGCCTTCGATACGGGA    EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAA/EEEEEEEEEAAAAA    NM:i:1  MD:Z:56G19  MC:Z:76M    AS:i:71 XS:i:76 RG:Z:../H2O XA:Z:chrM,-9379,76M,0;

```

Based on mapq and positions (respect to 569930) I see that the read with `cigar=18S58M` is skipped by mpileup.

**WHY?** The soft clip is way far from position chr1:569930

# Answer
> 3 votes
Examining the mpileup output for a single read on its own is a good way to figure out why bases are not appearing as expected.

So if you suspect that the `cigar=18S58M` read is the missing one (you are correct), prepare a SAM file containing just that read:

```
@SQ SN:chr1 LN:248956422
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE<A/EE//AE/AA///AAE///EEEEE//6///E//EE/6EE/EEEEEEEA/A/A/E/E//AEAEEEEE6A6AAA    NM:i:0  MD:Z:58 MC:Z:59M17S AS:i:58 XS:i:53RG:Z:../H2O  XA:Z:chrM,-9340,18S58M,1;

```

and run it through `mpileup`:

```
$ samtools mpileup -Q0 -B oneread.sam

```

This will show you that this read's base quality dips below 20 at the position you are interested in. Thus, as finswimmer suspected, your `-Q20` is removing it at this position.

# Answer
> 3 votes
Why do you believe it's the read with `cigar=18S58M`? The mapping qualities given in sam files are different to those in the mpileup output. I guess samtools is doing some recalculation.

You are using the parameters `-q20` `-Q20`. So you skip all reads with mapping quality \<20. This will result in 7 reads got to mpileup. Furthermore you are skipping bases with base quality \<20. Have a look at the base qualities for the read starting at position `569904`. If I count it correct you end up with a quality of `/` at position `569930` which means `14`. So I guess this read is dropped.

BTW: `samtools mpileup` is deprecated (only when creating bcf/vcf)  and should be replaced by `bcftools mpileup`.

---
Tags: samtools, mpileup
---


# Question
Title: Cluster is split in 2-3 locations on tsne plot - Suerat
I am running a single cell dataset (count data - exon) through Seurat. After running tsne I see a cluster (13) split in 3 different locations on the plot. Here are the commands I am running:

```
library('Seurat')
## rows (genes) that have 99% zeros were filtered out,
## and then columns (cells) that have 99% zeros were filtered out
isp.data <- read.table(file = "genes.cells99.percent.filtered", sep =
            ",", row.names = 1, header=TRUE)
isp <- CreateSeuratObject(raw.data = isp.data, project =
       "brain_regions")
isp <- NormalizeData(object = isp, normalization.method =
       "LogNormalize",scale.factor = 10000)
isp <- FindVariableGenes(object = isp, mean.function = ExpMean,
       dispersion.function = LogVMR,do.plot = FALSE)
hv.genes <- head(rownames(isp@hvg.info), 1000)
isp <- ScaleData(object = isp, genes.use = hv.genes, num.cores = 3)
isp <- RunPCA(object = isp, pc.genes = hv.genes, pcs.compute = 100,
       do.print = TRUE, pcs.print = 1:5, genes.print = 5)
PCElbowPlot(object = isp, num.pc = 100)
isp <- JackStraw(object = isp, num.replicate = 100, display.progress =
                 FALSE, num.pc = 30)
JackStrawPlot(object = isp, PCs = 1:30)
isp <- FindClusters(object = isp, reduction.type = "pca", dims.use =
        1:26, resolution= 0.75, print.output = 0, save.SNN = TRUE,
        force.recalc=TRUE)
isp <- RunTSNE(isp, dims.use = 1:26, do.fast = T)
TSNEPlot(isp, do.label = T, pt.size = 0.5)

```

Here is the tsne plot:

I have tried varying parameters for resolution in FindClusters (from 0.5-3.5), but that just increases the clusters found. Further I tried to play around with parameters min.dist (0.75,1), n\_neighbors (10 to 30) and perplexity(5 to 50) but I still get at least 1-2 clusters split in the tsne plot.

Another thing I have tried is to reduce the dimension (dims.use=10) on which tsne is ran but the clusters seem to be too close together, like so:

Note: Dimensionality of dataset (after filtering) - 23824 genes across 8967 samples.

My question being why does such behavior happen? Is my filtering strategy not correct? Have I not normalized/scaled data correctly?

# Answer
> 3 votes
Your cluster labels come from graph clustering implemented in the `FindClusters()` function. The resulting clusters are then visualised with a 2D tSNE plot (via `RunTSNE()` and `TSNEPlot()`). So, your cluster 13 is not split into three sub-clusters but cells within cluster 13 ***look*** somewhat ***distant*** from each other on the tSNE plot. Having mentioned distances on a tSNE plot, please see this SE Data Science entry Can closer points be considered more similar in T-SNE visualization? and especially the comment ***"What I mean is that you can't just use distance in the lower space as a similarity criterion"***. In your case "the lower space" would be your tSNE plot.

As far as I can see from your code, you are following the widely-adopted "Seurat workflow" and all seem fine including your TSNE plot output. I could not spot the QC step though: apart from the one at object initiation, one usually applies filters for number of genes and % mitochondrial gene contribution (using `FilterCells()` in Seurat v2 and `subset()` in Seurat v3).

It would be a good idea to try UMAP instead of tSNE, please see Evaluation of UMAP as an alternative to t-SNE for single-cell data (Becht et al, 2019). Seurat v3 vignettes are now using UMAP and not tSNE any more.

Lastly, please see How to Use t-SNE Effectively to learn about points to take into account when using/interpreting tSNE plots. You will see that you would get surprisingly dissimilar tSNE plots from the same underlying data with different parameters (some are hardcoded in the `RunTSNE()` function).

---
Tags: rna-seq, seurat, single-cell, clustering, normalization
---