# Question
Title: Folded Protein Chunk Dimensional Classification?
Are there known dimensional measurements for the classification of folded proteins given a starting chunk/domain as defined by something like the clustering functionality of MSM Builder?

Examples of what I would be looking for would be dimensions such as:

* pH
* Salinity
* Dissolved Sugars
* Reducing Agent Concentration
* Ambient Electronegativity

With corresponding scales and precisions such as:

* logorithmic vs linear
* step size
* distribution (e.g. dissolved/flat vs fractal/inclusions)

The above parameters would be included in folding software after a known starting conformation provided by MSM Builder or similar chunking software to derive a GUID for proteins without (or with extremely minimal and tuneable) collisions?

For instance, if hundreds of GUIDs on average end up equating to the same folded protein that's fine, but if one GUID could map to two or more thermodynamically stable folded proteins at greater than some incredibly small (e.g. 0.0000001%) then it wouldn't fit what I'm looking for. I know this is inherently NP-hard, but the complete ID would be something like:

```
{
    sequence,
    chunkId,
    pH,
    salinity,
    dissolvedSugar,
    reducingAgentConcentration
}

```

So I would expect the chunkId and dimensional parameters to play a significant role in reducing collisions for different starting parameters.

It's worth noting I'm not seeking some magic bullet where you type an ID and get a fully folded protein, just a reliable system of GUIDs with enough starting information embedded in the ID to start the folding process.

# Answer
> 3 votes
## CATH

The CATH database classifies protein by fold: https://www.cathdb.info/ So the value from that is probably the most useful for you.

## Crystallographic conditions

* pH
* Salinity
* Dissolved Sugars
* Reducing Agent Concentration

Your crystallising conditions do not mean much. They are a solution that is close to precipitating your protein but slowly enough for it to crystallise. A protein can crystallise under multiple different conditions and a crystal may be chosen simply because it is a lower index well in the Hampton screening conditions... These mean very very little. And unless you are planning on doing ML to determine the best screening condition these should be avoided.

## Data from the PDB

You have probably noticed that the only stats PDB/PDBe give are crystallographic, like space group, unit cell dimensions, resolution, completeness in shells, Rfree, Rmerge etc. So nothing there is of much use... except for bound ligand. Often this is the cofactor, but often it's a cofactor analog or a product (as opposed to the substrate). But that information can come from Uniprot etc.

## pH

A solution has a pH, a protein has a pI. At a pH equal to the pI it will be neutral —and less insoluble. ProtParam can calculate the pI from a linear sequence. However, surface charge totally shows different spots. APBS can map the solution to the Poisson-Boltzmann equation to your protein. There is a PyMOL wrapper, a Pyrosetta wrapper and more. But you still would need to some stats and 3D operations to convert it into a small vector of values —rototranslations and projects are not for the faint of heart or in a rush.

## Melting temperature

BRENDA has lots of enzymatic data on enzymes including melting temperature or kinetics at different temperatures (which follow the Eyring equation (TST), until MMRT kicks in) which could be useful.

For a given protein you can calculate the *relative* folding ∆∆G (eg. with Rosetta Score), which is a proxy to melting temperature in a buffer it is happy in.

Despite paining me to say... some fanciful Poisson-Boltzmann map derived vector and ∆∆G are probably all you could get from analysing structures. Whereas the rest are from CATH classification and other databases that are not structure based.

---
Tags: proteins, protein-structure, identifiers, thermodynamics
---


# Question
Title: How to find the number of contigs produced? N50 length in base pairs? N90 length in base pairs?
I am having trouble with some underlying questions about my project.

I have ran a Trinity tool to create a trinity.fasta file.

To determine some underlying questions, I used the utility asm\_stats to gather some basic statistics.

I am just having some issues of calculating data.

DATA:

Main genome scaffold N/L50: 7464/1.032 KB
Main genome contig N/L50: 7464/1.032 KB
Main genome scaffold N/L90: 14314/106
Main genome contig N/L90: 14314/106

Number of contigs (scaffold length):

25,534 (ALL)
25,534 (100)
21,609 (250)
10,818 (500)
907 (1 KB)
38 (2.5KB)
5 (5 KB)
2 (10KB)

With the data attached above, would my N50 and N90 length in **BASE PAIRS** so excluding the kb, be 7464 and 14314, respectively?

*Would the number of contigs produced by this assembly* be 7, since there are 7 different contigs produced or would it be 25,534?

Thanks!

# Answer
> With the data attached above, would my N50 and N90 length in BASE PAIRS so excluding the kb, be 7464 and 14314, respectively?

Looks like this is inverse so your N50 is 1032bp and your N90 is 106bp. N90 cannot be longer than N50 by definition. You L50 and L90 are 7464 and 14314

> Would the number of contigs produced by this assembly be 7, since there are 7 different contigs produced or would it be 25,534?

25534

Now for the real answer. **Do not calculate or report N50 on transcript assemblies. It is more than useless.** Transcriptomes are - by nature - dynamic in length. Longer is not better in this regard. Are your short transcripts lower quality assemblies, probably not, but N50 treats them as such. Please read the Trinity wiki for other, better ways to QC your transcript assembly.
> 1 votes

---
Tags: rna-seq, fasta, fastq, file-formats
---


title: kraken2 to OTU table creation

text:
# Question
Title: kraken2 to OTU table creation
I did run kraken2 and I get two kinds of output How do I generate OTU table from these outputs for phyloseq usage?

# Answer
> 1 votes
You need to run Bracken to the Kraken2 report output to estimate abundance.

---
Tags: genome, microbial-genomics
---


title: Contamination on genome assembly

text:
# Question
Title: Contamination on genome assembly
I had a question for the community.

I have a genome of a new species that has been sequenced via 150pb Illumina paired-end.

To verify the quality of the assembly I used the `BUSCO pipeline` which allows me to search for 99% of the BUSCO genes (gene not duplicated) present in a taxon in the genome of my species.

Globally I only miss `18%` of the BUSCO genes, which means that the assembly is not too bad.

Now when I plot the **G+C%** and the **coverage** of all the scaffolds containing BUSCO genes (so a priori non-contaminating) I should find a fairly homogeneous cloud.

Instead I end up with two different coverage profiles (noted red and blue below):

I then thought that the low coverage sequencing profile around 4x might be a eukaryotic contaminant like a acarians, however when I blast the BUSCOs present on these low coverage scaffolds I get hits that are very close to the genome of my species...

Does anyone have any idea what these are?

here is an output of the `contigsTable.csv` where the two first row are scaffold with low coverage in the figure, and two last are scaffold with high coverage in the figure

FRC output

Coverage distribution

Here is an Icarus summary of the assembly :

```
Assembly
# contigs (>= 0 bp) 1038631
# contigs (>= 1000 bp) 120132
# contigs (>= 5000 bp) 22245
# contigs (>= 10000 bp) 9986
# contigs (>= 25000 bp) 1322
# contigs (>= 50000 bp) 97
Total length (>= 0 bp) 655530043
Total length (>= 1000 bp) 451497861
Total length (>= 5000 bp) 260894112
Total length (>= 10000 bp) 174183054
Total length (>= 25000 bp) 45006012
Total length (>= 50000 bp) 5956145
# contigs 214542
Largest contig 120620
Total length 519138286
GC (%) 32.94
N50 5062
N75 1624
L50 21982
L75 71030
# N's per 100 kbp 22.28

```

Here is the BUSCO summary as well :

```
2544 Complete BUSCOs (C)
2312 Complete and single-copy BUSCOs (S)
232 Complete and duplicated BUSCOs (D)
1060 Fragmented BUSCOs (F)
811 Missing BUSCOs (M)
4415 Total BUSCO groups searched
C:57.7%[S:52.4%,D:5.3%],F:24.0%,M:18.3%,n:4415

```

EDIT

I ran the fastq sequences with another aligner (MEGAHIT) to see if it was not because of an assembly issue and I got the same profils coverage see here :

# Answer
> 5 votes
This could be organism-specific. We don't have a lot of info so far, so I would check a few more things:

Run something like FRC\_align. Check if there's a clear signal between regions flagged as suspicious by it and your coverage graph.

Is it a eukaryote? Plant? Check where mitchondria and chloroplasts are on the plot. They will have different GC/coverage signals than the rest of your assembly, which is fine, expected even.

Worried about contamination? Run kraken2 on your raw data and see what you get there.

Run RepeatMasker and check if repetive regions contribute to the effect you see.

Finally, unless your specimen is low heterozygosity, haploid, or inbred, you probably have some effect from ploidy. There are numerous ways to check that. A self-self dot plot with MUMmer and looking coverage histograms are a good start. The idea is that you should see two peaks (let's assume diploid). One for the split haplotigs, and one for collapsed.

I checked the GC distribution of 45 genomes in your order. It's kind of all over the place, but there are some bimodal ones:

What does GC vs Length look like? Here's the one from the 45 genomes mentioned above:

# Answer
> 0 votes
That is indeed puzzling. My best guess would be contamination, even if they carry some of the BUSCOs. I run all de-novo genomes through blobtools, which creates a similar plot to the one you made manually, but also adds taxonomical annotations. It can look like this:

Also, what is the proportion of the genome that is low-coverage? Does it have a significant span? I almost always see a bunch of tiny contigs with low coverage in de novo assemblies, but they usually represent a very small fraction of the genome.

-- edit --

One more thought, are you sure you don't have a tetraploid species? The two coverage peaks kind of correspond to 8 and 32x right? One is 1/4 of the other. We developed some tools to investigate this kind of problems in data.

-- edit 2 --

Just a few more comments. Coverage differences can have various reasons, but all the sequencing biases would cause a skew of the distribution, hence the apparent bimodality of your coverage (more apparent from the two 2d plots) suggests two distinct sources. In general, could be due to various ploidy levels or contamination either by bacteria, endosymbionts or even other species getting into DNA isolation. Now that polyploidy is ruled out, it boils down to contamination, and to figure out what contamination it is, I recommend the tool I liked above.

---
Tags: assembly, read-mapping, coverage
---


# Question
Title: about two group comparison of three group data in DESeq2 package
When we have three groups samples (A,B,C) with 25000 genes and the main interest is A vs B, should we limit samples only A and B for normalization to perform DEG analysis? Or better to include all samples (A,B, and C) for normalization and specify the comparison group within three groups?

```
res <- results(dds, contrast=c("tissu_type","A","B"))

```

# Answer
> 1 votes
As explained on the comments, having some other samples might distort your estimations, add an outlier or something alike.

The advantage of having all the three (or all) groups of samples before testing differences, because usually methods like limma, DESeq2 and others uses the information of variation of the genes across all samples to update their prior probabilities and the variation of the genes on the statistical model used to compare between groups.

If the **method** and your **model** are robust use all the three groups (generally the methods used are robust, but sometimes the models are not). By robust I mean that it doesn't miss any know (and available) source of variation. Otherwise the comparison might be comparing only the samples not the effect of the condition on the samples. If you only compare the samples it would be equivalent to use a t.test to the limma methods (not completely sure on the DESeq2 case). Or explained differently, you wouldn't take advantage of the statistical improvements implemented on these methods

Usually the reasons why these groups are sequenced on the same batch is because they are related (or to add some known ground truth), so you can release it with the data when communicating your results but you may omit that you used this third group on the methods section.

# Answer
> 0 votes
Plot out the PCA of your samples. If, say, PC1 splits the C's away from the others, and PC1 is more than, say, 50% of the variance, or if the C samples are a lot more spread out than the A and B samples, keep it apart. Otherwise, do them all together.

If the C's aren't too different from A and B, the math will benefit from including them.

---
Tags: r, rna-seq
---


# Question
Title: removing nextera transposase adapters, cutadapt
I am learning about NGS analysis and im currently learning about QCing and removing adaptors.

I am working on SRR1972920\_1.fastq file.

When running fastqc tool on that file, adapter contamination is present in the form of Nextera Transposase adapters.

Searching for that adapter sequence via google the sequence for at the 3' end to be removed is TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG (obtained via illumina adapter sequence-Nature).

Running a very basic cutadapt command line to remove the adapter

```
cutadapt -a TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG -o <any output name> <my input file>

```

After running the fastqc tool again on the new output file from cutadapt, it seems nothing changed and the contamination still present.

However, when running the same command above but changing the sequence to CTGTCTCTTATA, which is the sequence of transposase adaptors according to fastqc tool adapters\_list.txt file all is good after running fastqc tool on the file

```
cutadapt -a CTGTCTCTTATA  -o <any output name> <my input file>

```

Now, why this is happening? which Nextera transposase sequence is the correct? or am I doing something wrong?

# Answer
> 5 votes
The adapter sequence you have googled is the sequence on the adapter primer. This works when you want to remove primer-dimers.

With transposase adapters or ATAC seq, you have very short fragments and you sequence into the adapter. So what you is sequenced is actually the reverse complement. See below image from trimmomatic manual

So if you look at what worked, and the sequence you googled, they are reverse complement of one another.

You can also see this in the NexteraPE-PE.fa, Trans1 and Trans1\_rc @Haci recommended:

```
more NexteraPE-PE.fa
>PrefixNX/1
AGATGTGTATAAGAGACAG
>PrefixNX/2
AGATGTGTATAAGAGACAG
>Trans1
TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
>Trans1_rc
CTGTCTCTTATACACATCTGACGCTGCCGACGA
>Trans2
GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
>Trans2_rc
CTGTCTCTTATACACATCTCCGAGCCCACGAGAC

```

To answer your question, if most of your fragments are short, and there's not a lot of primer dimers, then yes you are doing the right thing.

# Answer
> 3 votes
I cannot quite help(\*) on the problem you are having with `cutadapt` but can point you out to `Trimmomatic`, for which the developers have been granted permission to distribute Illumina adapter sequences.

If you download `Trimmomatic`, you will see a bunch of Illumina adapter sequence files in the `/trimmomatic-0.39/adapters/` folder. Actually the adapter sequence you have found via googling is listed in the file `NexteraPE-PE.fa`.

* Regarding the undesired output with `cutadapt`, you were probably removing one of the multiple Nextra adapter options but not the exact one used in your library preparation step. You won't have to worry about finding the exact adapter sequence used in the library prep step when using `Trimmomatic` as it scans for the presence of all possible adapter sequences.

# Answer
> 2 votes
I came across the same problem of Nextera transposase contamination in my shotgun metagenome sequence. I specified the library in trimmomatic and Nextera transposase adapters were successfully removed. Please see the code below:

```
java -jar /opt/software/Trimmomatic/0.39-Java-1.8/trimmomatic-0.39.jar PE -phred33 1004_R1.fastq.gz 1004_R2.fastq.gz ../QC_data/1004_R1.fastq.gz ../QC_data/1004.qcup_R1.fastq.gz ../QC_data/1004_R2.fastq.gz ../QC_data/1004.qcup_R2.fastq.gz ILLUMINACLIP:/opt/software/Trimmomatic/0.39-Java-1.8/adapters/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

```

The output showed completed successfully with the following message:

```
TrimmomaticPE: Started with arguments:
 -phred33 1004_R1.fastq.gz 1004_R2.fastq.gz ../QC_data/1004_R1.fastq.gz ../QC_data/1004.qcup_R1.fastq.gz ../QC_data/1004_R2.fastq.gz ../QC_data/1004.qcup_R2.fastq.gz ILLUMINACLIP:/opt/software/Trimmomatic/0.39-Java-1.8/adapters/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
Using PrefixPair: 'AGATGTGTATAAGAGACAG' and 'AGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTCCGAGCCCACGAGAC'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTGACGCTGCCGACGA'
ILLUMINACLIP: Using 1 prefix pairs, 4 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 18205918 Both Surviving: 10467310 (57.49%) Forward Only Surviving: 7470543 (41.03%) Reverse Only Surviving: 48503 (0.27%) Dropped: 219562 (1.21%)
TrimmomaticPE: Completed successfully

```

And the QC report showed the Nextera transposase was successfully removed. \[!\[QC report\]\[2\]\]\[2\]

---
Tags: adapter, cutadapt
---


# Question
Title: "perl: warning: Setting locale failed." in RepeatMasker
I'm trying to run Repeatmasker in Linux on the command line with:

```
./RepeatMasker -species human -alu -gff -dir /mnt/lustre/users/Analysis/RepeatMaskerOutput /mnt/lustre/users/Software/STAR/bin/Linux_x86_64/HRG_Newest_and_Annotation/GCA_000001405.28_GRCh38.p13_genomic.fna

```

but I am getting this error:

```
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_GB.utf8",
        LC_ALL = "",
        LC_CTYPE = "en_GB.utf8",
        LANG = "en_GB.utf8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
RepeatMasker version open-4.0.9
TRF program not configured! This version of RepeatMasker requires a
local installatoin of TRF.  Please visit:
               http://tandem.bu.edu/trf/trf.html
to obtain the current version.  Once installed please re-run the
RepeatMasker configure script to setup RepeatMasker to use the new
installation. at ./RepeatMasker line 442.

```

I have downloaded rmBlast and the latest TRF file. I have put the TRF file in my `Users/username/bin` PATH folder and used `chmod a+x` (and `chmod 755`) on it. I then changed my `LD_LIBRARY_PATH` using `export` to my glibc directory.

After this I used `perl ./configure` and directed the program to the PATH folder where the TRF file was. It gave a series of the same error messages as above but was 'successful' as shown here:

```

[user@nodec23() RepeatMasker]$ perl ./configure
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
RepeatMasker Configuration Program

**TRF PROGRAM**

  This is the full path to the TRF program.
 This is now used by RepeatMasker to mask simple repeats.

Enter path  [ /mnt/lustre/users/Software/RepeatMasker/RepeatMasker/TRF.pm ]

Add a Search Engine:
   1. Crossmatch - Phrap             [ Configured ]
   2. RMBlast - Improved NCBI Blast  [ Configured, Default ]
   3. HMMER3.1 & DFAM                [ Un-configured ]
   4. ABBlast/WUBlast                [ Un-configured ]
   5. Done

Enter Selection: 5
 -- Setting perl interpreter...
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").

Congratulations!  RepeatMasker is now ready to use.
The program is installed with a the following repeat libraries:
  Dfam database version Dfam_3.0
  RepeatMasker Combined Database: Dfam_3.0
Further documentation on the program may be found here:
  /mnt/lustre/users/Software/RepeatMasker/RepeatMasker/repeatmasker.help

```

I'm at a loss of what to do next. I have googled and tried changing each of the locale settings with export to the following:

* en\_GB
* en\_US.utf8
* en\_GB:2000
* en\_GB.utf8

but none have changed the error.

I am running a Windows 10 machine and using PuTTY to connect to an HPC which uses a Linux environment.

Locale settings are as follows:

```
locale: Cannot set LC_CTYPE to default locale: No such file or directory
locale: Cannot set LC_MESSAGES to default locale: No such file or directory
locale: Cannot set LC_ALL to default locale: No such file or directory
LANG=en_GB
LC_CTYPE="en_GB"
LC_NUMERIC="en_GB"
LC_TIME="en_GB"
LC_COLLATE="en_GB"
LC_MONETARY="en_GB"
LC_MESSAGES="en_GB"
LC_PAPER="en_GB"
LC_NAME="en_GB"
LC_ADDRESS="en_GB"
LC_TELEPHONE="en_GB"
LC_MEASUREMENT="en_GB"
LC_IDENTIFICATION="en_GB"
LC_ALL=

```

Also please note that I do not have administrator access.

Typing 'perl -v' gives the following:

```
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").

This is perl, v5.10.1 (*) built for x86_64-linux-thread-multi

Copyright 1987-2009, Larry Wall

Perl may be copied only under the terms of either the Artistic License or the
GNU General Public License, which may be found in the Perl 5 source kit.

Complete documentation for Perl, including FAQ lists, should be found on
this system using "man perl" or "perldoc perl".  If you have access to the
Internet, point your browser at http://www.perl.org/, the Perl Home Page.

```

Typing "LC\_ALL=C perl ./configure" and going throught the configurations leads to the same outcome but without the perl errors:

```
Enter Selection: 5
 -- Setting perl interpreter...

Congratulations!  RepeatMasker is now ready to use.
The program is installed with a the following repeat libraries:
  Dfam database version Dfam_3.0
  RepeatMasker Combined Database: Dfam_3.0
Further documentation on the program may be found here:
  /mnt/lustre/users/k1894998/Software/RepeatMasker/RepeatMasker/repeatmasker.help

```

However, the error remains when trying to run Repeatmasker.

# Answer
> 0 votes
I think it is quite simple you need to install `perlbrew` here

You then run the software under perlbrew it should work.

At present you are running the system perl5 and that is causing all sorts issues, for example the age of the installation and the fact system operations and your software are using the same version of perl.

All the instructions are in the website link.

---
Tags: rna-seq, perl, linux, repeatmasker
---


# Question
Title: Is there an R package that computes homoplasy excess ratios (HER)?
I want to calculate the homoplasy excess ratio for a phylogenetic tree, for example homoplasy resulting from recombination.

I am wondering if there are any R packages out there that will do this?

I have looked at phyloseq, Rphylopars, phylocurve, PhylogeneticEM, PhyloinformR, but these packages do not have functions that calculate HER. I searched "homoplasy excess ratio" on bioconductor, but it returned ca. 3000 packages.

# Answer
> 1 votes
Maynard's Smith's homoplasy is what you want and it is available in the recombination detection package (RDP) described here. Its really a Windows package rather than R.

---
Tags: r, phylogenetics, phylogeny, covid-19, sars-cov-2
---


# Question
Title: Extract sequences from MultiFASTA aligned file, by coordinates
I am trying to extract a specific sequence from a multifasta file, from each sequence in the aligned file.

The sequences look like this, and there are 32 sequences within the multiFASTA:

```
>NC_045512.2_Severe_acute_respiratory_syndrome_coronavirus_2_isolate_Wuhan-Hu-1_complete_genome
ATT AAA GGT TTA TAC CTT CCC AGG TAA CAA ACC AAC CAA CTT TCG ATC TCT TGT AGA TCT GTT CTC TAA ACG AAC TTT
AAA ATC TGT GTG GCT GTC ACT CGG CTG CAT GCT TAG TGC ACT CAC GCA GTA TAA TTA ATA ACT AAT TAC TGT CGT TGA
CAG GAC ACG AGT AAC TCG TCT ATC TTC TGC AGG CTG CTT ACG GTT TCG TCC GTG TTG CAG CCG ATC ATC AGC ACA TCT
AGG TTT CGT CCG GGT GTG ACC GAA AGG TAA GAT GGA GAG CCT TGT CCC TGG TTT CAA CGA GAA AAC ACA CGT CCA ACT

```

Using the following code I can pull the sequence, using coordinates, from a single FASTA file, but how would one iterate this to pull the section (by coordinates) from each sequence within the multiFASTA?

```
from Bio import SeqIO
from Bio.SeqIO.FastaIO import FastaWriter
import os

os.chdir("C:/Users/")

record = SeqIO.read("Whole genome.fas", "fasta")

with open("Spike protein.fasta", "wt") as out:
    writer = FastaWriter(out, wrap = [])
    writer.write_header()
    writer.write_record(record[21563:25384])
    writer.write_footer()

```

Thanks in advance!

# Answer
I developed an answer as follows:

```
from Bio import SeqIO
from Bio.SeqIO.FastaIO import FastaWriter
import os

os.chdir("C:/Users/")

with open("Spike protein.fasta", "wt") as out:
    for seq_record in SeqIO.parse("Whole genome.fas", "fasta"):
        out.write(str(seq_record.id) + "\n")
        out.write(str(seq_record.seq[21563:25384]) + "\n")

```
> 1 votes

# Answer
You can use `SeqIO` for the writing in addition to reading the sequences. This has the advantage that it will always write the correct fasta header, which is not always just the `.id` if there is whitespace, e.g. (untested)

```
from Bio import SeqIO

with open("Spike protein.fasta", "w") as out:
    for seq_record in SeqIO.parse("Whole genome.fas", "fasta"):
        seq_record.seq = seq_record.seq[21563:25384]
        SeqIO.write(seq_record, out, "fasta")

```

BTW don't put spaces in your filenames, it can lead to issues
> 2 votes

---
Tags: fasta, python, covid-19, sars-cov-2, seqio
---


# Question
Title: Does module size decrease as the number of samples increase in WGCNA?
We are working with a lot of samples, more than 50,000 people! We are using WGCNA to get modules and it seems like there are less number of module genes the more people we use. We are following the tutorial guidelines. We set the min modulesize = 30 and deepSplit = 4 (max sensitivity) However, we are getting very small modules with numbers almost 30 for all. We want to know if using this many samples creates an artifact and the variation between the samples may reduce the correlation seen between the genes? If that is the case, do you suggest any way we can fix that? Also do you think z-scaling expression data improves the module detection?

Thank you so much.

# Answer
> 1 votes
The biggest difference a large number of samples makes is that you can usually decrease the soft thresholding power. If you use power 6 (or 12 for a signed network), try decreasing it to 3 or even 2 (6 or 4 for signed network). Apart from that, it could also be that the co-expression modules in your data are simply smaller than what people usually find, nothing apriori wrong with that.

---
Tags: wgcna
---


# Question
Title: Is it possible to convert BAM file from one genome assembly to the other?
I Have multiple BAM files that are referenced to UCSC genome assembly GRCh37/hg19 that are read in different time frames. Now, I am planning a different studies that require assembling all the data against one assembly. Will it be possible to convert all the BAM files that I have to an upgraded version or convert the assembly to a different new assembly and how? I am new in this and not a bioinformatician so, if you have any questions regarding my question to answer you can ask, any type of help is highly appreciated.

# Answer
It is possible, with tools like CrossMap, but my recommendation is to re-align. Assemblies can change a lot between versions and your reads will undoubdedly align differently. Re-aligning is cheap, it will take a few days and you're good. Re-doing your analysis because you realize the coordinate translation messed something up will cost you much more time and energy.
> 1 votes

---
Tags: sam, genome, assembly
---


# Question
Title: Extracting concatenated metadata from a column in R
I've downloaded some ChIP-seq data from www.chip-atlas.org and I want to use this to conduct some enrichment tests with genomic regions from my experiment. The .bed file generated by that website is set up so that each row represents a genomic region from a specific ChIP-seq experiment. I have the data set up as a dataframe in R.

Each row has a single column with all of the metadata, including such important information as the target transcription factor, model system, and tissue type. These are delimited by semi-colons, but the number of pieces of information varies between rows so it's impossible to simply split by the delimiter.

The metadata is structured something like this:

```
name=CTCF; strain=B6; age=adult; tissue=brain
name=Mecp2; age=embryo; tissue=5t3 cells; passages=5-10

```

So within each row, the order and number of metadata terms varies. I want to reshape this column into something like this:

```
name; strain; age; tissue; passages
CTCF; B6; adult; brain; NA
Mecp2; NA; embryo; 5t3 cells; 5-10

```

How can I extract and reshape this data in R? It seems like it should be possible, but the normal approaches to splitting columns by delimiters are failing me since the rows all have different numbers of items and order.

Here's a slice of the real data which is much more messy than my little example.

```
test_df <- data.frame(pos=c(3002983, 3002881, 3010946, 3021775), metadata=c("ID=SRX661585;Name=Ctcf%20(@%20Brain);Title=GSM1446329:%20CTCF%20brain%20ChIPSeq%3B%20Mus%20musculus%20x%20Mus%20spretus%3B%20ChIP-Seq;Cell%20group=Neural;<br>source_name=Mouse%20adult%20hybrid%20brain;tissue=Adult%20whole%20brain;strain=BL6/spretus;chip%20antibody=CTCF,%20MillIpore%2007-729;", "ID=SRX2957352;Name=Smc1a%20(@%20AtT-20);Title=GSM2684855:%20SMC1%20ChIPseq%20in%20AtT-20%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;<br>source_name=AtT-20%20cells;cell%20line=AtT-20%20cells;cell%20type=pituitary%20corticotroph%20cell;passages=10-20;chip%20antibody=Bethyl%20A300-055A;", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;<br>source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;<br>source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);"))

```

Thank you for your help!

# Answer
In general you want to:

1. Clean your text
2. Split the features by semicolons
3. Split keys and values by '='
4. Create a temporary data frame with keys and value
5. Reshape

```
library(tidyverse)
library(urltools)
library(textclean)
test_df <- data.frame(pos=c(3002983, 3002881, 3010946, 3021775), metadata=c("ID=SRX661585;Name=Ctcf%20(@%20Brain);Title=GSM1446329:%20CTCF%20brain%20ChIPSeq%3B%20Mus%20musculus%20x%20Mus%20spretus%3B%20ChIP-Seq;Cell%20group=Neural;<br>source_name=Mouse%20adult%20hybrid%20brain;tissue=Adult%20whole%20brain;strain=BL6/spretus;chip%20antibody=CTCF,%20MillIpore%2007-729;", "ID=SRX2957352;Name=Smc1a%20(@%20AtT-20);Title=GSM2684855:%20SMC1%20ChIPseq%20in%20AtT-20%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;<br>source_name=AtT-20%20cells;cell%20line=AtT-20%20cells;cell%20type=pituitary%20corticotroph%20cell;passages=10-20;chip%20antibody=Bethyl%20A300-055A;", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;<br>source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;<br>source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);"))

tbl <- apply(test_df, 1, function(row){
  # Split by smicolons and clean html tags
  n <- replace_html(unlist(strsplit(row[2], ";")))
  # Decode urlencoded characters, mostly space
  n <- url_decode(n)
  # Split into keys and values
  s <- do.call(rbind, str_split(n, "="))
  # Get keys and values into vectors and clean up leading and trailing whitespace
  v <- s[, 2]
  v <- str_replace_all(v, "(^ +)|( +$)", "")
  k <- s[, 1]
  k <- str_replace_all(k, "(^ +)|( +$)", "")
  # Create final data frame with original columns (pos), and also columns from
  # keys and values
  tibble(pos = row[1], key = k, value = v)
})

# Concatenate list into tibble
(tbl <- do.call(rbind, tbl))

```

This tibble looks like:

```
# A tibble: 37 x 3
   pos     key           value
   <chr>   <chr>         <chr>
 1 3002983 ID            SRX661585
 2 3002983 Name          Ctcf (@ Brain)
 3 3002983 Title         GSM1446329: CTCF brain ChIPSeq; Mus musculus x Mus spretus; ChIP-Seq
 4 3002983 Cell group    Neural
 5 3002983 source_name   Mouse adult hybrid brain
 6 3002983 tissue        Adult whole brain
 7 3002983 strain        BL6/spretus
 8 3002983 chip antibody CTCF, MillIpore 07-729
 9 3002881 ID            SRX2957352
10 3002881 Name          Smc1a (@ AtT-20)
# … with 27 more rows

```

Using this, you can:

```
# Pivot wider
tbl %>%
  pivot_wider(names_from = key, values_from = value)

```

Which nets you:

```
# A tibble: 4 x 15
  pos   ID    Name  Title `Cell group` source_name tissue strain `chip antibody` `cell line` `cell type`
  <chr> <chr> <chr> <chr> <chr>        <chr>       <chr>  <chr>  <chr>           <chr>       <chr>
1 3002… SRX6… Ctcf… GSM1… Neural       Mouse adul… Adult… BL6/s… CTCF, MillIpor… NA          NA
2 3002… SRX2… Smc1… GSM2… Neural       AtT-20 cel… NA     NA     Bethyl A300-05… AtT-20 cel… pituitary …
3 3010… SRX5… Mecp… GSM3… Neural       Olfactory … NA     NA     MeCP2 antibody… NA          Olfactory …
4 3021… SRX5… Mecp… GSM3… Neural       Olfactory … NA     NA     MeCP2 antibody… NA          Olfactory …
# … with 4 more variables: passages <chr>, `strain/background` <chr>, `genotype/variation` <chr>, age <chr>

```
> 1 votes

---
Tags: r, chip-seq, data-preprocessing
---


title: Multiple Alignment cost application

text:
# Question
Title: Multiple Alignment cost application
Is there any biological application/case where someone would be interested in the estimated total cost of the alignment between a set of sequences (genes or amino acids) without the aligned sequence itself?

Edit: to clarify by cost I mean the cost function such as the edit distance. I have some ideas for heuristics on the multiple sequences alignment problem (which is intractable) but these would not retrieve the aligned sequence, therefore, I was wondering whether such an algorithm would be useful.

# Answer
**UPDATE**

In response to comments, it seems that the question is more along the lines of:

**"what are some applications where computing summary statistics associated with a specific MSA result could be useful, even if we can't return the sequences?"**

I have 2 responses to this:

1. in keeping with my original answer, I would say that it would still be useful to obtain some measure of identity, similarity, or pairwise distances between sequences (or any summary statistic about the relationship of the sequences). There are lots of applications for this level of data. I don't know whether this is something that can also be extracted from the computation, I'm a little unclear on what exactly the proposed computation is doing if it's getting a cost but not sequences.
2. Michael says below (and I agree) that an alignment cost is likely to strongly depend on particularities of the underlying alignment model as well as the data, and things get nasty when you try to account for gaps or weirdness. Nonetheless, it's possible that a cost could be used as a sort of guide heuristic for a separate follow-up alignment procedure, to indicate a portion of parameter space which is likely to yield good results. I don't know what this looks like, but it could be simple as resource provisioning for the computation (e.g. CPUs/RAM/disk). Alternately, maybe the cost result can be a little richer and actually seed some parameters in the next computation.

Again, without more detail I'm not sure what this looks like (and I'm not an alignment expert so my suggestions wouldn't be very good anyways).

3. The simplest application might be to come up with a quick "check" of the sequences to ask whether they are alignable in the first place. Like, as a QC step (did you forget to reverse complement?) maybe. (this assumes that the alignment in question is pretty onerous and you don't want to have to try several times.)

In summary, a lot depends on how you get that cost and what else you *could* get in addition to the cost. I wonder if rather than focusing on cost you can output some other richer summary of the alignment (such as e.g. sequence similarities, or gapping stats, etc.).

**ORIGINAL ANSWER**

The short answer is "yes".

Quite frequently we are interested in measuring "identity" or "similarity" of two sequence strings. As one totally arbitrary example, in metagenomics, quite frequently two rDNA amplicons that show 97% identity are considered to be in the same species, 95% same genus, 98% same strain, etc. So these measures are frequently useful as heuristics.

Identity itself is just a length-normalized edit distance across a sequence, after all.

But more generally, such measures of identity or similarity can themselves be understood as cost measures, in the sense that e.g. a dynamic programming alignment algorithm is trying to find an optimal set of matches between two sequences, which necessarily will maximize identity. Generally what it is actually maximizing is similarity (according to some cost matrix such as e.g. BLOSUM or Jukes-Cantor or whatever).
> 1 votes

# Answer
I'm not a fan personally, but it is context specific. The problem is that it is heading towards optimisation, particularly algorithmic optimisation. If the alignment is the final calculation - that is cool, but the situations are select. It depends how far you want to take it.

The rationale of my skepticism is simple - if an alignment was decided on a 'maximization' criteria, be it likelihood or some other metric - the subsequent calculations based on the alignment, such as phylogenetics Jukes Cantor, Kimura 2 and 3-parameter etc ... we enter circular mathematics. The critiera that defined the alignment is the criteria we are using to perform unrelated post hoc analysis. In truth there is a degree of circular mathematics in any algorithm used to produce an alignment, but we overcome this by optical verification precisely avoid this.

In this sense the alignment is simply to establish homologous amino acid/nucleotide positions for downstream analysis. Thus optical optimisation, whilst cautioned against in every other area of bioinformatics, is considered useful in alignments because of the risk that a logical error in the alignment is amplified in downstream analysis, potentially leading to an erroneous conclusion.

Algorithms are a problem where indels are concerned and aligning across an indel. It is very easy to make errors here and 'optical analysis' would simple introduce a blockwise deletion to avoid the uncertainty, i.e. if you are not sure just chop it out - better to lose data than risk artefacts on an equivocal judgement call.

I'm not saying don't try, I am saying don't look to overstate the applicability of the metric.
> 2 votes

---
Tags: phylogenetics, phylogeny, gene, sequence-analysis, multiple-sequence-alignment
---


# Question
Title: How do you currently track bugs in python using freely available software/websites?
I am starting a project using Python to track, archive, assign and manage bioinformatics bugs. This software will be able to be integrated into project management software or stand alone.

* How do you currently track bugs/assign bioinformatics bugs using freely available bioinformatics software/websites?
* What are features that you like/dislike? What are features that you would like to see in a bioinformatics bug tracker?

I am a Biology student with basic/core Python skills.

# Answer
> 2 votes
Use IPython and activate the trace in the Python debugger as follows

```
%pdb 1

```

or

```
%pdb on

```

to turn this off `%pdb 0` or `%pdf off`

I do recommend using a formal IDE such as PyCharm or Visual Studio Code. The professional version of PyCharm is highly rated and within that use **breakpoints** . In Python3.7 there is supposed to be `breakpoint()` so you may not need an IDE, i.e. you set the breakpoint in your code using `breakpoint()`. In older versions of Python you can import this function via, `import pdb; pdb.set_trace()`

Finally check out `PDB++` which is found here

Its important to note that experienced programmers often debug by simply inserting `print()` throughout their code and don't bother with an IDE

---
Tags: python, bugs
---


# Question
Title: Is there a simple command for outputting a tab delimited columns?
I am working on a fasta file and am writing my command in nano within command-line and executing using python, also within a command line.

My objective is to get my command to provide me with a tab delimited file with three columns: first column should contain my **sequence name**, second column should provide me with my **sequence length**, and the third column should show the **sequence itself**.

I have written the following command so far within nano:

For example, I would like my command to provide me with the desired output and with the following order: Gene name ; Gene length ; Gene seq

```
A06842X 45 GORHQRIHQEERHEUWOHRPPTRWFAWWEAKJNFWEJQEFQEPRT

Jb4329PC 21 LDHKQAJTYOWEPLKFREEQW

```

# Answer
> 1 votes
I have no idea what do you want to do with fasta files reformatted like this, but I think something like this should do the job (untested):

```
from Bio import SeqIO
import sys

for olig_fasta in SeqIO.parse(sys.argv[1], "fasta"):

  name = olig_fasta.name
  seq = str(olig_fasta.seq)
  seq_len = str(len(olig_fasta.seq))

  print("\t".join([name, seq_len, seq]))

```

---
Tags: python, fasta, assembly, file-formats
---


title: What type of Protein ID is this?

text:
# Question
Title: What type of Protein ID is this?
I got a FASTA file from NCBI RefSeq. These proteins are annotated with IDs like this : protein\_id=XP\_019970915.1. What type of ID is this? I would like to have them converted in TrEMBL IDs in order to do, later, a GO term enrichment but I have to know where I am starting from. Thank you in advance

# Answer
I'm not sure I understand the question, but XP\_019970915.1 is an 'NCBI Reference Sequence' ID for 'erythrocyte binding antigen-181' (https://www.ncbi.nlm.nih.gov/protein/XP\_019970915.1/). If you look that protein up in UniProt you can find other identifiers/info and GO annotations etc (https://www.uniprot.org/uniprot/A0A151LVZ9). Does that answer your question?
> 2 votes

---
Tags: go-enrichment, identifiers
---


# Question
Title: How can I subset WGS data to the level of WES variants?
I would like to compare mutational signatures1 in patients from different studies, however some studies are based on exome seq (i.e. ~20,000 coding variants) and some are from whole genome seq (i.e. ~22,000 coding variants) - is there a way to 'downsample' WGS data to better reflect WES data and effectively 'ignore' the coordinates of those ~2000 coding variants in the VCF files?

> 1Alexandrov, L.B., Kim, J., Haradhvala, N.J., Huang, M.N., Ng, A.W.T., Wu, Y., Boot, A., Covington, K.R., Gordenin, D.A., Bergstrom, E.N. and Islam, S.A., 2020. The repertoire of mutational signatures in human cancer. Nature, 578(7793), pp.94-101.

# Answer
> 3 votes
bcftools would be my choice, I'm sure bedtools could do the trick, too. Something along this

```
bcftools view --regions-file

```

or --targets, --regions-file might require a tabix index.

# Answer
> 2 votes
I'll expand slightly on the previous answer. First print off the positions for the exome file and then use bcftools view to filter the variants from the whole genome file. You could also index the whole\_genome.vcf file to make the filtering faster.

```
bcftools query \
    -f'%CHROM\t%POS\n' \
    exome.vcf > exome_variants.txt

bcftools view \
    -T exome_variants.txt \
    whole_genome.vcf > whole_genome.exome_positions.vcf

```

---
Tags: genome, ngs, vcf, software-recommendation, exome
---


# Question
Title: Is there a Python program that is able to show convergent and divergent evolution on a phylogenetic tree?
I was wondering if there was a Python program that detects convergent and divergent evolution on a tree. I am also curious if it is just better to look at a phylogenetic tree and conclude what is convergent or divergent evolution than what a program tells you.

# Answer
> 1 votes
Well in the absence of a definition of 'divergent evolution' I would simply use a definition of allopatry. If the divergence is congruent with geographic distance or geographic foci this provides some proof of divergent evolution.

In this context you are looking for something like the biogeographic BAyesian analysis such as BayArea or geographic add-ons to `Beast`

There isn't a singular Python library that will handle this but there is some functionality with Biopython here

```
from Bio import Phylo

```

---
Tags: python, phylogenetics, gene, database, genetics
---


# Question
Title: Merging several fq.gz files or R1 and R2 classes into a single one
I have the following files:

```
V300066187_L4_B5RDBATtnuRAAAAA-408_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-408_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-408_2.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-408_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-407_2.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-407_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-405_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-406_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-405_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-407_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-407_1.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-406_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-406_1.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-405_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-405_1.fq.gz

```

the ending defines whether it is forward \_1 or reverse \_2.

I want to merge all the `_1.fq.gz` files into like `Merge1_.fq.gz` and the same for `_2.fq.gz`

I tried doing `cat *_(number).fq.gz > merged(number).fq.gz` and when I used the merges for assembly with MaSurca, it didn't detect any reads and failed. I used one pair and it worked well, but the genome size was smaller than expected, real size = 2 GB, estimated = 1.8GB.

Would appreciate how to solve this issue and use all my short reads for assembly - I need to input one file for all forward reads and one file for all reverse reads.

# Answer
> 2 votes
What you did should work. You don't need to decompress and recompress.

I'd check the files. Were they trimmed? Can you get the top, say 1000 lines of the merged files to work? (You will need to decompress to count lines accurately)

---
Tags: genome, fastq, bash
---


# Question
Title: How to remove low abundance and less prevalent data from my dataset?
I have a relative abundance table in .tsv format where samples are in column and rows contain the features (pathways). Something like this reproducible example:

```
# Pathway   ERR2601 ERR26014    ERR260145_Abu   ERR260  ERR2601 ERR2601 ERR2607 ERR2600 ERR26
1CMET2-PWY: N10-formyl-tetrahydrofolate biosynthesis    0.00826631  0.00946472  0.00962571  0.00411153  0.00811879  0.00844421  0.00915602  0.00789428  0.0091798
3-HYDROXYPHENYLACETATE-DEGRADATION-PWY: 4-hydroxyphenylacetate degradation  0   0   0   0.000695608 0   0   0   0   0
7ALPHADEHYDROX-PWY: cholate degradation (bacteria, anaerobic)   0   0   0   0   0   0   0   0   0
AEROBACTINSYN-PWY: aerobactin biosynthesis  0   0   0   0   0   0   0   0   0
ALLANTOINDEG-PWY: superpathway of allantoin degradation in yeast    0.0000331927    0   0.0000335389    0   0   0   0   0.0000189167    0.0000908934
ANAEROFRUCAT-PWY: homolactic fermentation   0.00682776  0.00681683  0.00675732  0.00495716  0.00691969  0.0053589   0.00514634  0.00779249  0.00426403
ANAGLYCOLYSIS-PWY: glycolysis III (from glucose)    0.00866693  0.0100437   0.00829714  0.00524566  0.00960262  0.0075093   0.00716347  0.00916205  0.00866289
ARG+POLYAMINE-SYN: superpathway of arginine and polyamine biosynthesis  0.00106395  0.000544088 0.000774645 0.00370237  0.00102368  0.00133848  0.00148691  0.00114062  0.00188739
ARGDEG-PWY: superpathway of L-arginine, putrescine, and 4-aminobutanoate degradation    0   0   0   0.00164625  0.0000278133    0   0   0   0
ARGININE-SYN4-PWY: L-ornithine de novo biosynthesis 0.00166327  0.00383993  0.000105528 0.00230283  0.00404469  0.00115975  0.000460711 0.00101816  0.00367102
ARGORNPROST-PWY: arginine, ornithine and proline interconversion    0.000132331 0.000244069 0.0000796126    0.000479006 0.000237704 0.000290783 0   0.00066458  0.0000581229
ARGSYN-PWY: L-arginine biosynthesis I (via L-ornithine) 0.0124159   0.00955758  0.0118644   0.00533077  0.00774959  0.0098352   0.0127232   0.00996349  0.0108839

```

Now, I want to keep features (i.e. pathways) which are present with abundance \>0.0001 and present in at least 10% of samples. What bash command can I use in order to do that?

# Answer
> 1 votes
Here is a bash/awk solution:

```
awk -F'\t' 'NR==1{print}; NR>1{N=0;for(i=2; i<(NF+1); i++){if($i >= 0.0001) N++}; if(N>0.1*NF) print}' input.tsv

```

We make sure the input field separator is tab, then we loop on the number of fields (NF) and add to the tally the `N` of fields above threshold. At the end we simply `print` the entire line if `N` is high enough.

EDITED: for the header, we print the first line (`NR==1`) as is, then apply the expression when Number of Records \> 1.

# Answer
> 2 votes
Assuming you have this stored in a matrix called `m`, where the pathway labels are row names:

```
IDX = apply(m, 1, function(x) sum(x > 0.0001)/length(x) >= 0.1)
m_subset = m[which(IDX),]

```

---
Tags: bash
---


title: RAD Seq Data Analysis without barcode

text:
# Question
Title: RAD Seq Data Analysis without barcode
I received a ddRAD Seq data from my supervisor without barcodes and restriction enzymes. I asked him for both and he said I don't need it since the data has been cleaned by the company. Now, I am kinda confused and don't know where to start.

# Answer
> 2 votes
Wow Biobash, that's kinda harsh. It is cool that sometimes you unexpectedly learned something new. What Michael mentioned is pretty interested, no one has ever mentioned that or I have never heard that.

So, my approach would be that you now are doing a Poolseq study, since you can't identify individuals. The other way to try to catch the paralogues is to map your reads to a sister species. My question is what are you trying to estimate? Genetic differentation? Genetic diversity? Pop Gen structure? https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6116966/ https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.5240 https://www.g3journal.org/content/9/12/4159

One semi-robust approach to assess usefulness is thinking about your system and defining what you expect to find, let's say in terms of genetic differentiation or genetic structure, and estimating that from your data. If that hypothesis is confirmed, your dataset might be useful. You can test the same by first estimating the loci by mapping your reads to a sister reference genome, to try to detect paralogues. Since you don't have barcodes, then each "locus" is independent, and you lose power, but there is nothing else you can do.

Good luck;

# Answer
> 0 votes
The fundamental problem with RAD is that if you hit paralogue - you're done.

Identifying a paraologue isn't necessarily easy depends on the species. Essentially if you're doing pop gen you need homologues and RAD has too many situations where it can fall over without you knowing.

What happens if you've a paralogue is that you start confusing concerted evolution (pop gen of a genome) with population structure between individuals. I do agree that without barcodes identifying individuals - thats a bit tricky ;-) , however RAD an approach loaded with pitfalls.

Pop gen means being able to identify populations and on face value that will not be trivial without barcodes (I suspect they are still there), but RAD is far from an ideal tool and can give weird and wonderful metrics which have little to do with the population structure.

Its not the only problem with RAD for the purpose you are using it.

---
Tags: snp, phylogenetics
---


title: Random GWAS data generator

text:
# Question
Title: Random GWAS data generator
I was wondering if there is a tool/script/program that randomly generates GWAS data. The purpose of such a tool would be to use it for educational purposes. So you generate some random .ped, .map and other filetypes. Various distributions like the MAF of all the SNPs, X-Chromosome classificaiton error rates and so on could be set before generating the data. General stuff one can then filter out with plink or similar software to increase the quality of the data.

Then a student would analyze these random files in plink and I already know what the outcome should be, therefore I can see right away if he or she did it right.

Does anyone know if such a program already exists? Or is there a somewhat easy way to program this myself?

Maybe it works if I use existing gene-data from NCBI where I use the given MAF to generate n variations the alleles based on that? Just some random SNPs from dbSNP?

# Answer
> 1 votes
Agreed with the above I would use an existing dataset from NCBIs dbSNP here

Random doesn't work well because HW signal will be meaningless as soon as the Jukes-Cantor distance is exceeded (which it will) and plink needs homologues to generate any signal either in 'reasonable' panmixia or clonality. The result will be panmictic beyond even the most panmictic organism ever known. I don't mind having a simulation of 100% clonalty, because 100% clonality does occur, either through inbreeding or binary division, but this simulation would generate the complete opposite result.

---
Tags: snp, gwas, plink
---


# Question
Title: Survival analysis using CoxPH - Effect of covariates
Hi and sorry for the long post in advance,

I'm doing a survival analysis of lung cancer patients using Python's lifelines package. According to the documentation, the function plot\_partial\_effects\_on\_outcome() plots the effect of a covariate on the observer's survival. I do understand that the CoxPH model assumes that the log-hazard of an individual is modelled by a linear function of their covariates, however, in some cases the effect of these covariates seems to perfect when plotted. I'll illustrate these with three images below describing the issues I have in my data analysis.

Example 1: This one I'm providing as it looks realistic to me, there is a clear difference in the effect of the covariate when positive or negative.

Example 2: This one just looks gimmicky, my concern is firstly there seems to be a perfect relation between the pathologic stage of cancer and the survival, it's logical that higher stages should impact the survival negatively, however, I don't understand why the graph is so symmetric as the survivability of cancer shouldn't be proportional to the patient cancer staging.

Example 3: Same problems as for example 2, this time with a binary not categorical variable (so I'm less reluctant on disbelieving it), however, the curves do seem to symmetric, almost like it's the same curve but with a standard deviation. I also highlighted a region of interest which I don't understand why both curves stabilize (maybe lack of patients in those regions?)

# Answer
> 3 votes
Example 2 and 3 seem symmetric because they are forced to by the model. The Cox proportional hazard relies on the proportional hazards assumption. This basically means that the ratio of the hazards are constant over time. If you were to plot the hazard ratio, it would look like a flat line over all `Event_at`. Another way of saying this is that the hazard ratios are constant over time.

How that shows up in a survival plot is a little different because of the relationship between hazards and survival. The survival curves end up following the pattern exemplified in Example 2.

As for example 3, I would guess that there are few individuals with event times past 3500. This is why you have those large jumps.

**Proportional Hazards**

So what does the proportional hazards assumption mean? Below are some plots of hazard functions to help your intuition on what the proportional hazards assumption means. In all of the following scenarios, the hazard ratio is 2 for blue-to-red. However, the proportional hazards assumption is only true in A, B, and C.

Since the *proportion* is remaining constant, it can look a little hard to see. To make it easier on the eyes, you can log-transform the hazard functions. Below is a plot that does that.

As you can see, the hazards are a constant distance from each other for A, B, and C. Basically, the proportional hazards assumption is saying that your scenario looks something like one of those plots.

In scenario D, the hazard ratio varies over time. At the beginning it is basically HR=1. At the end it is basically HR=4. With the proportional hazards assumption, Cox regression will essentially take the 'average' of the hazards over time. In this case it would give you an estimate of HR=2. Depending on the scenario, this 'smoothing' may be more or less acceptable. But this is basically what the proportional hazards assumption does in the background.

The proportional hazards assumption is actually more problematic in the cases where the two hazard functions cross each other. Below is an example of that.

In this scenario, blue has an increased hazard at the start (relative to red) then a decreased at the end. A Cox proportional hazards model would give you HR=1 in this case, because the average over all time points is no relationship.

So how do you tell if you violate the proportional hazards assumption? You may want to use a stratified Kaplan-Meier by your variables to get an idea of what the hazards look like without stipulating the proportional hazards assumption. You can use the Kaplan-Meier estimator to get a non-parametric estimate of the discrete hazard (i.e. it does not make the proportional hazards assumption). Then you can log-transform those hazards to get a plot like the second. There are also statistical tests and other diagnostics (Schoenfeld residuals) to check this assumption.

The next question is do those violations matter? This somewhat depends on the question you are trying to answer and how substantial the violations are. In cases where the two log-hazard curves don't cross (like in D), it may be acceptable to smooth over the variations. However, it settings where the curves cross in the final figure, it is probably best to not ignore this problem.

The original write-up of this answer was here

---
Tags: python, visualization, cancer
---


# Question
Title: If a gene is expressed at a level of 1/1200 compared to the average gene, how is probability 50:50 that we have a read mapped to it?
I am reading a book about RNAseq analysis and it says

> "To calculate the probability that a read will map to a specific gene, we can assume an average gene size of 4000 nt (100 M nt divided by 25,000 genes). At 30 M reads equivalent to 30× coverage, at single read 100 nt (or paired-end read 50 nt) length, we can expect a single read to map to the average expressed and length gene 4000 nt× 30 coverage/100 nt 1200 times. Thus, if the gene is expressed at a level of 1/1200 compared to the average gene, then we have a 50:50 probability to have a read map to it."

I understand that a single read would map to the average expressed 1200 times but I don't understand how the 50:50 probability is calculated. Please could you explain how this probability is determined.

# Answer
I actually disagree with that.. I guess I write this down as a discussion.

The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read.

However, because of sampling, sometimes you get 0, 1 or 2. If we assume the sampling follows a poisson distribution, we can calculate the probability of getting 0, 1 , 2 ... reads:

```
plot(dpois(0:10,lambda=1),ylab="Probability",xlab="No of reads")

```

The probabilities are like:

```
data.frame(n=0:5,p=dpois(0:5,lambda=1))

  n           p
1 0 0.367879441
2 1 0.367879441
3 2 0.183939721
4 3 0.061313240
5 4 0.015328310
6 5 0.003065662

```

The probability of getting exactly 0 or 1 read should be ~0.37. The probability of getting at least 1 read would be 1 - P(n=0) = 1- 0.367879441 = 0.6321206
> 3 votes

---
Tags: rna-seq, reads, coverage
---


# Question
Title: Measure the purifying selection for certain taxa along a phylogeny
I am posting this message because I need clarity about the analysis I want to perform. In my analysis, I have a homologous gene in **13 species**, and I would like to evaluate the selection pressure of this gene (its dN/dS) in **3 species** of the tree and test if this value is significantly different from the neutral model (dN/dS = 1).

So here is the protocol I am considering, and I would like you to tell me if I am right or not:

To calculate this selection pressure, the ratio of non-synonymous to synonymous substitution rates (dN/dS) must be calculated.

For this I use CODEML implemented in ete3.

So I need :

```
A phylogenetic tree of the gene (performed with the AA alignment of the gene via Iqtree)

A codon alignment (performed with the DNA alignment of the gene with pal2nal)

Mark branches in the tree.

```

Here is the exemple where I want to get the **dN/dS** value of **SP5,6** and **7**.:

```
         /-SP1
        /-|
       |   \-SP2
       |
     /-|      /-SP3
    |  |   /-|
    |  |  |   \-SP4
    |   \-|
    |     |   /-SP5
  /-|      \- |  /-SP6
 |  |         \-|
 |  |           |
 |  |            \-SP7
 |  |
 |  |   /-SP8
-|   \-|
 |      \-SP9
 |
 |     /-SP10
 |  /-|
 | |   \-SP11
 |-|
   |   /-SP12
    \-|
       \-SP13

```

And I want to calculate the **dN/dS** of `SP5,SP6` and `SP7` in order to confirme the hypotesis that the gene is under strong purifying selection on this 3 species.

I mark the SP5, SP6 and SP7, run a `free-branch` model and look at the Codeml output (the dN/dS value should be the same for the 3 marked species right?) Then this value is kind of the mean dN/dS value between three leafs?

Suppose I get an omega = **0.15** If I want to test that this value is significantly different from 1, I need to mark the same sheets and run the neutral model, then compare the two models (`neutral` vs. `free`) using chi2, and if the p-value is `<0.5`, then the dN/dS value (**0.15**) is significantly different from 1. I'm right ?

# Answer
> 1 votes
The following is not the clearest explanation of dN/dS humankind has seen, but I hope you get the gist and core idea of likelihood is absolutely correct. I am wondering if you are thinking of a different test to the way I would do this calculation. This ain't a text book clarity answer.

Anyway what you are proposing is a likelihood ratio test (LRT) for selection based on the parameterisation of omega. Your idea isn't quite correct ... its not dN/dS value you are looking at to assess significance, its the likelihood. This point is critical and is the most important difference in our thinking on this test.

LRT works because its a log likelihood and subtracting two log likelihood is the same as dividing them (I forget the theory - but you subtract the likelihood to work out the p-value via Chi squared). There's some theory about why the Chi-squared distribution works but I can't remember it.

You need to work out the degrees of freedom which is the parameter differences between 'free branch model' and the neutral model and this value is LOW. If the neutral model is the "global model", then its 2 degrees of freedom - IN YOUR TREE. This is because the SP you are looking at form a single clade (monophyly). From memory you stick a #1 on the branch leading into the clade, ie. defining the node and thats your 'model' on the treefile. Here you have two dN/dS one for the whole tree excluding your clade and one for your taxa. This gives you a likelihood.

**Note** you could stick a #1 on each branch within the clade, i.e. the treefile. That means something different - like every branch was subject to ... in this case purifying selection ... not just the common ancestor.

You now work out the likelihood of the 'global' model... i.e. remove the #1 and make sure it gives you only 1 dN/dS value.

Subtract the likelihoods and multiple by two - if its LRT (super important).

You now need to work out the degrees of freedom. If its the same model but one has an additional node with a different dN/dS ... the df = 2. If you've dumped a #1 on every branch in the clade ... I still think df = 2 - you need to check.

If the difference in likelihood results in *p* \< 0.05 you are correct the difference is significant ... whatever that difference might be.

If the neutral model is some different calculation of omega within the population genetics process (see here) then it will be the parameter differences between your model and the 'neutral' model.

**To reiterate** I like simulations so much I type in reiterations ... So again the LRT usual will usually subtract one likelihood from another and THEM multiply the probability by two ... and compare the result on the Chi squared table for the difference in degrees of freedom. And again what you are looking at is not the value of dN/dS but the difference in LIKELIHOOD. That is very important. You need to be careful with the degrees of freedom on an LRT because they can introduce an additional subtraction constant, test dependings (its ages since I've done one so I can't remember).

If the test is significant and one dN/dS is x for the clade and y for the rest of the tree, then that differences is significant and if y is \>1 (bit unlikely) and x \< 1 then sure thats purifying selection. Usually clades evolve via positive selection rather than purifying selection.

From memory that is my understanding ... but if you have a neutral model of omega and a purifying model of omega ... yeah sort of, but usually thats for a 'global' model, i.e. no differences within the tree. In your case your interested in a specific clade. The idea is still the same its measured through likelihood.

Its so long since I've done this, but I do caution not to confuse two different tests.

---

I do understand Yang et al did an MCMC for dN/dS ... what I'm thinking of here is maximum likelihood. You'd need to watch for that carefully and thats what you are referring to when you say Chi-squared, because that is LRT 100%. LRT is maximum likelihood.

---
Tags: sequence-alignment, phylogeny, positive-selection
---


# Question
Title: Does BLASTClust guarantee that proteins in different clusters are dissimilar?
I need to find dissimilar proteins. Looking through the PDB I found the weekly BLASTClust results of proteins that are 30% *similar*. However, I do not know if protein A in cluster 1 is guaranteed to be less than 30% similar to protein B in cluster 2.

The best and only documentation I can find for BLASTClust is here (link).

It states:

> The program begins with pairwise matches and places a sequence in a cluster if the sequence matches at least one sequence already in the cluster.

That's good, but consider a set of three proteins (A, B, and C). A is 40% similar to B and 40% similar to C. B is 20% similar to C. Say that the algorithm places B in cluster 1 and then C in cluster 2 as they are dissimilar. What happens when it tries to place A? Because A is similar to both B and C will they be joined into one cluster (this would be the best result for me) or is A just placed into one cluster even though it is similar to a protein in another cluster?

# Answer
> 2 votes
Based on the documentation (that I linked to in comment), it seems that the clustering algorithm used by blastclust is single-linkage hierarchical clustering.

SO what you can expect is that A will be placed with the protein to which it is most similar. And then the last protein will be added as an outgroup to that cluster due to its similarity to A (assuming that there are no other proteins).

So if your goal is "I want clusters that are guaranteed to contain dissimilar proteins", then you might not be able to naively take blastclust output and interpret it as such. You could either postprocess it to find the set of clusters in which no protein is more than 30% similar to a protein in any other cluster (which may just be one big cluster!), or find another clustering algorithm that is more suited to your task.

**Update**

It does look like blastclust uses a minimum *similarity* (not identity) threshold (from blastclust docs). So they use an implementation of single-linkage clustering that stops before merging everything into a single cluster:

```
If the coverage is above a certain threshold
 AND
the score density is above a certain threshold,

these two sequences are considered to be neighbored.

Thus determined neighbor relationships is considered symmetric and provides
the base for clustering by a single-linkage method (which puts a sequence
to a cluster if the sequence is a neighbor to at least one sequence in the
cluster)."

```

I don't think that this is a general property of the generic clustering algorithm though (see the "Naive algorithm" section).

---
Tags: blast
---


# Question
Title: How to remove batch effect from TCGA and GTEx data
I want to play some differentially expressed genes for gliomas cancer. I collected cancer tissue profiles from the TCGA database and normal tissue profiles from the GTEx database. I straightforward use ComBat from R package "sva", and the result shows little differentially expressed genes (DEGs) were found. How can I remove batch effects correctly?

Thanks

# Answer
> 1 votes
There is a paper that did what you are asking: Unifying cancer and normal RNA sequencing data from different sources.

# Answer
> 0 votes
This has been asked and answered here and Biostars. My take home message from many similar posts is "not to bother with correcting batch effects when performing differential expression analysis but to include 'batch' as a covariate". The latter accounts for the batch effects.

EDIT: This does not really answer OP's question, as @swbarnes2 and @ATpoint pointed out, they have a completely confounded design and the linear regression model would not accept such a model matrix.

---
Tags: combat, gtex
---


title: Chip-seq down-sampling problem

text:
# Question
Title: Chip-seq down-sampling problem
**Control sample chip read stat**

```
56234445 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
55173337 + 0 mapped (98.11% : N/A)
56234445 + 0 paired in sequencing
28078024 + 0 read1
28156421 + 0 read2
53808040 + 0 properly paired (95.69% : N/A)
54626508 + 0 with itself and mate mapped
546829 + 0 singletons (0.97% : N/A)
162718 + 0 with mate mapped to a different chr
112602 + 0 with mate mapped to a different chr (mapQ>=5)

```

**IgG sample chip seq read stat**

```
32745255 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
31546173 + 0 mapped (96.34% : N/A)
32745255 + 0 paired in sequencing
16362391 + 0 read1
16382864 + 0 read2
30934194 + 0 properly paired (94.47% : N/A)
31269484 + 0 with itself and mate mapped
276689 + 0 singletons (0.84% : N/A)
79786 + 0 with mate mapped to a different chr
49658 + 0 with mate mapped to a different chr (mapQ>=5)

```

**knockout sample Rep 1 stat**

```
28044731 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
27459222 + 0 mapped (97.91% : N/A)
28044731 + 0 paired in sequencing
14014483 + 0 read1
14030248 + 0 read2
26685254 + 0 properly paired (95.15% : N/A)
27162846 + 0 with itself and mate mapped
296376 + 0 singletons (1.06% : N/A)
116684 + 0 with mate mapped to a different chr
88965 + 0 with mate mapped to a different chr (mapQ>=5)

```

**Knockout sample rep 2 stat**

```
32410094 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
31774629 + 0 mapped (98.04% : N/A)
32410094 + 0 paired in sequencing
16196814 + 0 read1
16213280 + 0 read2
30802808 + 0 properly paired (95.04% : N/A)
31561218 + 0 with itself and mate mapped
213411 + 0 singletons (0.66% : N/A)
89688 + 0 with mate mapped to a different chr
63838 + 0 with mate mapped to a different chr (mapQ>=5)

```

Do i need to down-sample with respect to the IgG for control sample which has lot more reads at the same time Knockout sample 1 has less read. So is it good to go for further analysis or i need to down sample them ?

Any suggestion is really appreciated

# Answer
No, downsampling in general has no use as it reduces power. Peak callers such as `macs2` will internally scale data to account for sequencing depth differences.

The IgG itself is only useful for peak calling as it somewhat corrects for the unspecificity of the antibody to precipitate any kind of protein rather than the protein of interest. I hear people say, and personally agree with this by experience with my own experiments, that of course IgG pulls down notably fewer DNA that the target antibody, therefore the library requires notably PCR to produce enough material for sequencing. And this then makes it questionable how representative IgG is as a control. I personally prefer chromatin input to account for the intrinsic sequence/GC bias of the underlying genomic region in the sequencing/mapping. To my knowledge there are currently few established tools to really take these controls into account for downstream analysis such as differential testing via interaction terms because the library composition is simply so different from the IPs making normalization difficult.

Your focus imho should be rather to properly normalize your samples to correct for differences is signal/noise ratio, see for example my post over at biostars on how one could do that beyond naive per-million scaling.
> 1 votes

---
Tags: chip-seq
---


# Question
Title: How to concatenate two fasta sequences by py
I have three backbone and i want to concatenate 70 sequence into these backbone. such like: fasta file 1:

```
Backbone1: ACTGCGTATGTCAGCATATAGACCGCCT
Backbone2: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTC
Backbone3: GTCATGCTACGTGCTACGTGCTCGTCGCGC

```

fasta file 2:

```
sequence_1:CGCTAGTGTAGCTGCTACTCCG

```

output sequence:

```
Backbone1+sequence_1: ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
Backbone2+sequence_1: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG
Backbone3+sequence_1: GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCGs

```

# Answer
> 1 votes
First of all, those are not fasta sequences. Fasta is a specific format, and requires a `>` followed by a sequence identifier on one line and the one or more lines of sequence.

Assuming you really have two files in the format you show, you can do:

```
$ awk -F: '{ gsub(/ /,""); if(NR==FNR){bbone[$2]; next} for(b in bbone){print b$2}}' file1 file2
GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCG
ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG

```

The `gsub(/ /,"");` just removes the extra space you have in `file1` after each `:`. If you also want to have the sequence names as you show, use:

```
$ awk -F: '{ gsub(/ /,""); if(NR==FNR){bbone[$2]="Backbone"NR; next} for(b in bbone){print bbone[b]"+"$1":", b$2}}' file1 file2
Backbone3+sequence_1: GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCG
Backbone1+sequence_1: ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
Backbone2+sequence_1: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG

```

If you really need to do this in python for some reason, please edit your question and show us your python script so we have something to work on and expand. I really don't want to write a whole python program for something so trivial.

---
Tags: python, sequence-analysis, phylogenetics
---


title: Scooping and GeneMatcher

text:
# Question
Title: Scooping and GeneMatcher
A colleague has uploaded about a dozen candidate disease genes to GeneMatcher, described here and here and located here, with excellent success. However last year and again recently he has encountered the case wherein a case was uploaded and a cohort formed, but a few months down the line one of the dropped out collaborators published an n=1 paper in a low tier journal, basically using GeneMatcher to check whether the disease link of the gene was legitimate, to assess the phenotypic range and see what test are worth doing.

* Is there documented precedence of this behaviour (preferably citable)? If not,
* does anyone have personal experience of this and what is the ratio of incidences over total submissions (I know that is hearsay, but it would be very helpful)?
* And also has anyone successfully taken action (say getting user barred)?

# Answer
> 1 votes
I am having trouble finding anything that directly addresses the issue of failed collaborations leading to scooping, but here are some examples of similar behavior:

As for punishment, that sounds hard. Scooping is sort of built into science; fighting over precedence is a big deal and has been for a long time. My solution would be to get rid of current academic publication and have everyone's open digital lab notebooks streaming results to publication feeds in real time; that way there's no ambiguity. Preprints and micropublications are sort of a poor man's version of this.

This is not a very complete answer but there may be enough there to help.

---
Tags: gene, mutations, publishing
---


# Question
Title: I have used TER to break the long bonds of a chain in my PDB
I’m now not sure what I need to alter in my PDB to get it to work in leap. I know breaking the bonds turns the formerly connected residues into terminal residues. It keeps saying that 3 of my atoms no longer have a type , how can I fix this ?

# Answer
Normally you have three options in a PDB:

* Add a `TER` record
* Add the  `OXT` atom, if your chain was cleaved or cleaved spontaneously (*i.e.* Asp-Pro sequence), the C-terminal would sill have a terminal oxygen. i.e. making it a terminal residue. If you have protons, you'll need to add also `3H` .
* Assign a different chain

If you are using forcefields as is your case, the second is the only choice. The `C`  atom and the `N`  atoms have a type that does not match because they have the wrong number of bonds —hence the type error.
> 2 votes

---
Tags: rna, pdb
---


title: Moving UMI from tag into read name

text:
# Question
Title: Moving UMI from tag into read name
I have a bam file with Unique Molecular Identifiers (UMIs) for each read present on the RX tag ( such as `RX:Z:TGAGAAGGG`), as expected by picard and fgbio tools.

However, several other tools (such as UMI-tools) expect the UMIs to be present at the end of the read name, such as (`NB501943:73:HK333BGXC:3:12402:20313:16437_TGAGAAGGG`).

Does anyone know of a tool to do this? I have an `awk` solution which I'll self-answer with, but I don't know whether that is robust.

# Answer
My current solution is (`awk` is not a tool I'm very familiar with, so this may not be good `awk` code, there is presumably a way to get it to print the header lines unmodified.):

```
samtools view -H temp.sam > temp2.sam
samtools view temp.sam |
awk '{OFS = "\t"} {for(i=1;i<=NF;i++)  if ($i ~ /^RX:Z:/) {$1=$1"_"$i; gsub("RX:Z:","",$1); print}}' >> temp2.sam

```

This seems to work when I test it, but I don't know if I'm missing anything.
> 1 votes

---
Tags: sam, umi
---


# Question
Title: What datasets are available out there for prediction based on DNA sequences?
I am looking for publicly available data for a genomics deep learning project. My goal is to compare different architectures to predict biological insights from DNA sequences.

I have heard about Janggu before, but I am looking for a few datasets which are like the MNIST equivalents in genomics. Kaggle doesn't seem to have many such datasets, and I would be happy if you share your insights about potential datasets that one could use.

## For example, it could be datasets that map 100-Mers to their ChIP-seq values or datasets that annotate a DNA sequence with the probability that a transcription factor would bind to it.

Choong and Lee lists a few datasets on their IEEE paper from 2017:

# Answer
Good question and the Janggu link is cool. Thanks!

Kaggle is/was heavily image recognition. There isn't an MNIST type data set in genomics, in case anyone is wondering this is a classic natural language processing (NLP) type analysis to identify handwritten single digit numbers.

There has been several FDA precision challenges, one being described here:

FDA precision challenges

This could provide something towards a 'standard'.

A key issue remains the accurate vectorisation of sequence and whether you think that is/was satisfactory in my personal opinion.

There is an old review article here
> 2 votes

---
Tags: phylogenetics, phylogeny, deep-learning, tensorflow, pytorch
---


# Question
Title: Plot information stored in dataframe within a tree (ggtree)
Hello to the entire Stackoverflow community! I'm writing to you because I'm currently building a phylogeny with ggtree and I have in parallel a table like this one: **df** :

```
Node_number Number
50   31
1    15
13   25
35   12
40   4
5    2

```

and a tree such as :

and the idea is that I woule like to incorporate the df$Number information in each Node in the tree and get somethhing like I put in the picture (expected part)

tree data :

```
structure(list(edge = structure(c(31L, 32L, 33L, 34L, 34L, 35L,
35L, 33L, 32L, 36L, 37L, 37L, 36L, 38L, 39L, 40L, 40L, 39L, 41L,
42L, 42L, 43L, 43L, 41L, 44L, 44L, 45L, 45L, 46L, 46L, 47L, 48L,
48L, 47L, 38L, 49L, 50L, 51L, 51L, 50L, 52L, 53L, 53L, 54L, 54L,
52L, 49L, 55L, 55L, 31L, 56L, 57L, 57L, 58L, 58L, 56L, 59L, 59L,
32L, 33L, 34L, 1L, 35L, 2L, 3L, 4L, 36L, 37L, 5L, 6L, 38L, 39L,
40L, 7L, 8L, 41L, 42L, 9L, 43L, 10L, 11L, 44L, 12L, 45L, 13L,
46L, 14L, 47L, 48L, 15L, 16L, 17L, 49L, 50L, 51L, 18L, 19L, 52L,
53L, 20L, 54L, 21L, 22L, 23L, 55L, 24L, 25L, 56L, 57L, 26L, 58L,
27L, 28L, 59L, 29L, 30L), .Dim = c(58L, 2L)), tip.label = c("t5",
"t9", "t25", "t14", "t30", "t26", "t2", "t10", "t12", "t15",
"t1", "t3", "t6", "t23", "t17", "t18", "t4", "t22", "t27", "t16",
"t29", "t7", "t21", "t20", "t13", "t11", "t19", "t24", "t8",
"t28"), edge.length = c(0.912875924259424, 0.293603372760117,
0.459065726259723, 0.332394674187526, 0.65087046707049, 0.258016780717298,
0.478545248275623, 0.766310670645908, 0.0842469143681228, 0.875321330036968,
0.339072937844321, 0.839440350187942, 0.34668348915875, 0.333774930797517,
0.476351245073602, 0.892198335845023, 0.864339470630512, 0.389989543473348,
0.777320698834956, 0.960617997217923, 0.434659484773874, 0.712514678714797,
0.399994368897751, 0.325352151878178, 0.757087148027495, 0.202692255144939,
0.711121222469956, 0.121691921027377, 0.245488513959572, 0.14330437942408,
0.239629415096715, 0.0589343772735447, 0.642288258532062, 0.876269212691113,
0.778914677444845, 0.79730882588774, 0.455274453619495, 0.410084082046524,
0.810870242770761, 0.604933290276676, 0.654723928077146, 0.353197271935642,
0.270260145887733, 0.99268406117335, 0.633493264438584, 0.213208135217428,
0.129372348077595, 0.478118034312502, 0.924074469832703, 0.59876096714288,
0.976170694921166, 0.731792511884123, 0.356726912083104, 0.431473690550774,
0.148211560677737, 0.0130775754805654, 0.715566066093743, 0.103184235747904
), Nnode = 29L), class = "phylo", order = "cladewise")

```

does someone have an idea ?

# Answer
You can accomplish this using the ggtree package available on Bioconductor.

First you will need to combine your tree with the data.

```
library(tidyverse)
library(ggtree)

ftree <- tree$edge %>%
   as_tibble() %>%
   mutate(Node_number = 1:n()) %>% # finds edge numbering
   right_join(data, by = "Node_number") %>% # find internal node associated with your node numbering
   left_join(fortify(tree, ladderize = FALSE), ., by = c("node" = "V2")) # merge tree and data

```

Note that this will produce a warning, but it does work.

Then plot the tree using ggtree and its `geom_nodepoint` layer.

```
ggtree(ftree, ladderize = FALSE) + # plots tree
   geom_tiplab() + # adds tip labels
   geom_nodepoint(aes(subset = !is.na(Number), size = Number, x = branch)) + # adds your data as points
   scale_size_continuous(range = c(4, 12)) # scales point sizes

```

You can adjust the size of the points with the `scale_size_continuous` function.

Data

```
data <- read_csv('Node_number,Number
50,31
1,15
13,25
35,12
40,4
5,2')

```

and `tree` is the same as you provided.

Ggtree also has a decent manual available at https://yulab-smu.top/treedata-book/.
> 4 votes

---
Tags: r, phylogenetics, phylogeny, ggplot2, ggtree
---


title: How do you calculate the top L/5 score?

text:
# Question
Title: How do you calculate the top L/5 score?
I am currently investigating how to predict protein structure and contact map predictions. There, I see things like top L/5 and top L/2 scores as a way to evaluate the contact map. What formulas do you use to evaluate these? Also, what is a meaningful evaluation method?

# Answer
> 2 votes
A contact map predictor will usually output the probability of a contact for every residue pair in a protein, often ommitting very close pairs (e.g. less than 6 residue separation).

To calculate the L/n score, rank the predicted contacts descending by probability and look at the first L/n contacts. The fraction of these that are actually contacts in the native structure is the L/n contact precision score.

Often, you would filter the list to look at only contacts of a certain sequence separation, e.g. long range contacts between residues of \>23 separation (definitions vary between papers). If you have a good predictor this score is high for L/10 and falls as you go to L/1, which is expected as proteins don't usually have that many long range contacts to find.

---
Tags: proteins, protein-structure
---


# Question
Title: How to translate amino acid sequences to Nucleotide sequences
I want to convert a list of fasta ( protein sequences) in a .text file into corresponding nucleotide sequences. A Google search gives me result of DNA to protein conversion but not vice versa. Also, I came across How do I find the nucleotide sequence of a protein using Biopython?, but this is what I am not looking for. Is there any possible way to do it using python.Moreover, I would like to solve it using python programming. I am sure there must be some way to do it rather than writing a code from scratch. Thanks!

# Answer
Although there is not a unique nucleotide sequence that translates to a given protein, one can list all the possible DNA sequences that do translate to that protein.

An online tool that does just that is Backtranambig, from EMBOSS. It produces a DNA sequence representing all the nucleotide sequences matching the input protein, using IUPAC ambiguity codes.
> 5 votes

# Answer
DNA Chisel (written in Python) can reverse translate a protein sequence:

```
import dnachisel
from dnachisel.biotools import reverse_translate

record = dnachisel.load_record("seq.fa")
reverse_translate(str(record.seq))
# GGTCATATTTTAAAAATGTTTCCT

```
> 2 votes

# Answer
You can't do this because there is redundancy in the genetic code and the *same* protein sequences can be encoded by *different* nucleotide sequences. There are 64 codons and ~20 amino acids, e.g. `GCT`, `GCC`, `GCA` and `GCG` all encode Alanine.
> 1 votes

---
Tags: proteins, python, phylogenetics, sequence
---


title: How to choose docking score cutoff?

text:
# Question
Title: How to choose docking score cutoff?
I made a virtual screening of a large database (~5M) and want to filter ligand structures for further more accurate screening (from Schrodinger HTVS to SP resolution). My next step is to define a cutoff docking score, but it seems arbitrary to choose one. Are there any established methods for filtering the ligands out?

# Answer
> 3 votes
It is utterly arbitrary and depends on your resources.

We can all agree your cutoff would be a lot lower during the ordering step than at followup *in silico* analyses, but even then if you were ordering the top *N* at \\$100-250 each, should you spend \\$1,000 or \\$10,000?

In the *in silico* step you have real time and CPU-time to factor in making this call even blurrier.

But there are some things to put those scores into context:

* The units are kcal/mol, i.e. these are predicted ∆∆G of binding. Here is an explanation of what does a kcal/mol mean for docking — 1 kcal/mol is a hydrogen bond or the average collision energy of water at 37°C.
* If you have a crystal structure with a known ligand, re-docking and cross-docking that would give you a benchmark of what you should expect ——keeping my mind that natural ligand are horrid and don't obey Lipinski's rule. The natural ligand will wander and to assess how good the docking went RMSD is used, but for this I mean looking at its GlideScore as if it were any other compound.
* Rescore with a different tool —Gold etc.

## Ligand efficiency

An important thing to remember is that ∆∆G is size dependent —a huge peptidomimetic will score more than a smaller molecule that hits the hotspots perfectly. Therefore, ligand efficiency (LE), ∆∆G over heavy atoms (reviewed in detail here) is a good metric to prevent huge molecules dominating spuriously. If your library has a variety of different sizes, it is best to sort by LE than ∆∆G alone. In fact, if you get a few small hits you can merge/link them or expand on them, whereas if you get nothing because you gambled on larger hits you have to start from square one. I would argue that this is true for repurposing screens and the low success of repurposing VLS despite the enthusiasm is due to the neglect of this.

---
Tags: protein-structure, docking
---


title: How to fix NameError?

text:
# Question
Title: How to fix NameError?
I am working on a project and am having issues with the following code that I have written in `nano`:

```
from Bio import SeqIO
import sys
import re

    fasta_file = (sys.argv[1])
    for myfile in SeqIO.parse(fasta_file, "fasta"):
      if len(myfile) > 250:
       gene_id = myfile.id
       mylist = re.match('H149xcV\_\w+\_\w+\_\w+',gene_id)
       print (">"+list.group(1))

```

This is the error I receive when I execute my command on `command-line`:

```
File "mpo.py", line 7, in <module>
    gene_id = myfile.id
NameError: name 'myfile' is not defined

```

I have a fasta file with the format

```
>H149xcV_Fge342_r3_h2_d1 len=210 path=[0:0-206]
ACTATACATGAGGAGAACATAGAACAAAAATGGGACCATAGATATATAACAATAGAAGATATAGAGAACACAATAGACAACTTATTAGGAAAGAGGTGTGTCGTCATGGAGCTGATGTTCGAGGATACTTTGCATGGTCATTCTTGGATAATTTTGAGTGGGCTATGGGATACACCAAGAGGTTTGGCATTGTTTATGTTGATTATAAGAACGGGC

>H149xcV_ytR1oP_r3_h2_d1 len=306 path=[0:0-207]
ATTAGAGTCTGAGAGAGTCTTGATTTGTCGTCGTCGAGAAATATAGGAGATCTGATTAGAGGAGAGAGCGGCCTAGGCGATGCGCGATATAGCGCTATATAGGCCTAGAGGAGAGTCTCTCTCTTTTAGAAGAGATAATATATATATATATATGGCTCTCCGGCGGGGCCGCGCGAGAGCTCGATCGATCGATATTAGCTGTACGATGCTAGCTAGCTTATATTCGATCGATTATAGCTTAGATCTCTCTCTAAAGGTCGATATCGCTTATGCGCGCGTATATCG

```

I wish to reformat my file so that it only provides me with the unique gene id's for genes with a length greater than 250 bp.

My desired output is:

```
>H149xcV_Fge342_r3_h2
>H149xcV_ytR1oP_r3_h2
>H149xcV_DPN78333_r3_h2
>H149xcV_AgV472_r3_h2
>H149xcV_DNP733_r3_h2

```

# Answer
You might consider changing the code to make the whitespace a little more consistent (as it was, it failed when i tried to run it).

Also, why are you using regular expressions? You have already extracted the `gene_id` using Biopython. The following works for me:

(updated to handle the underscore delimiting)

```
from Bio import SeqIO
import sys
import re

fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile.seq) > 250:  # note .seq here!
        gene_id = myfile.id
        # unclear why using regex
        #match_list = re.match('H149xcV\_\w+\_\w+\_\w+',gene_id)
        #print (">"+match_list.group(0))
        print("_".join(gene_id.split("_")[0:-1]))

```
> 1 votes

# Answer
You wrote,

```
mylist = re.match(H149xcV\_\w+\_\w+\_\w+)

```

try...

```
myregex = re.compile(r'H149xcV\_\w+\_\w+\_\w+')
mylist = myregex.match(myfile.id)

```

or

```
mylist = re.match(r'H149xcV\_\w+\_\w+\_\w+',myfile.id)

```

Watch your indentation it should be 4 spaces. Note you are losing the last '\_\w+', but don't think it matters.

I think it is better to use,

```
mylist, _, _ = re.split(r' ', myfile.id)

```

The problem with your code is `list` will just return `true` ... and that is not what you want. If you are using `split` carefully watch the number of ' ' in your code you can also use `r'\s'`

You can also set maxsplit as a parameter, which might be useful.

---

`myfile` not defined.

Please try,

```
fasta_file = (sys.argv[0])

```

Generally your argv command isn't cool and not just \[0\] versus \[1\]. I forget the argv module everyone uses .. there's two of them, they are a bit fiddle but are very cool of options. You way in is okay, okayish, sort of okay ... Anyway you appear to be requesting the second input into the commandline, wherea you are probably placing the fasta file in the first position.

Also

```
print (fasta_file)

```

Check something happens that you picked up the file.

The other way in is,

```
mypath = '/User/username/location/fasta.fa'
for myfile in SeqIO.parse(mypath, "fasta"):

```

where `mypath` is a OS X style, you're probably using Linux so its... /home/username/pathtofile/fastafile.fa (something like that) or

```
for myfile in SeqIO.parse('/User/username/location/fasta.fa', "fasta"):

```

I forgot .... don't ever use `list` as a variable, by `mylist` is fine. Again I've not debugged my code.
> 3 votes

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