# Question
Title: GROMACS Coordinates and velocity and force trajectories do not match accuracy
I am using GROMACS to perform MD simulations of proteins. So, I tried to record the simulation results `nstxout`, `nstvout`, and `nstfout` with the same step size (`nstxout 1, nstvout 1, nstfout 1`). Then, I calculated the coordinate (step n) of an arbitrary step of the trajectory (step n) from the velocity of step n - 1, step n, and calculated the velocity of n - 1/2, and used that velocity to calculate the coordinate (n-1) of step n using the velocity of n - 1/2. LeapFrog was used in the simulation.

How can I get the coordinates, velocity, and force to be consistent?

# Answer
> 4 votes
From GROMACS documentation about reproducibility:

> The following factors affect the reproducibility of a simulation, and thus its output:
>
> * Precision (mixed / double) with double giving “better” reproducibility. Number of cores, due to different order in which forces are accumulated. For instance (a+b)+c is not necessarily binary identical to a+(b+c) in floating-point arithmetic.
> * Type of processors. Even within the same processor family there can be slight differences.
> * Optimization level when compiling.
> * Optimizations at run time: e.g. the FFTW library that is typically used for fast Fourier transforms determines at startup which version of their algorithms is fastest, and uses that for the remainder of the calculations. Since the speed estimate is not deterministic, the results may vary from run to run.
> * Dynamic linking to different versions of shared libraries (e.g. for FFTs)
> * Dynamic load balancing, since particles are redistributed to processors based on elapsed wallclock time, which will lead to (a+b)+c != a+(b+c) issues as above
> * Number of PME-only ranks (for parallel PME simulations)
> * MPI reductions typically do not guarantee the order of the operations, and so the absence of associativity for floating-point arithmetic means the result of a reduction depends on the order actually chosen
> * On GPUs, the reduction of e.g. non-bonded forces has a non-deterministic summation order, so any fast implementation is non-reprodudible by design.

I would add that in your case of hard coding everything by hand you may have to look a bit into your barostat and thermostat (although their influence should be found on speed of the atoms). Maybe a sample code with coordinates and speed of 4 atoms would help to give a more precise answer.

You can also try having a look at gromacs source code. I believe OpenMM provides a bit more control and clarity on each step, if for some unfathomable reason you are not doing this only for curiosity sake but because you actually have to make a code that reproduces these results.

---
Tags: simulated-data, gromacs
---


# Question
Title: How to draw a 3d structure using 3d coordinates from pdb file?
I am trying to create a 3d structure from 3d coordinated from a pdb file.

I am trying to draw 3d coordinates from pdb file "ATOM 1 N VAL A 14 22.692 -27.591 -39.402" and using webgl javascript.

I've managed to get the 3d coordinates from a pdb file but i don't know the exact logic that is used to create 3d structure from the coordinates found in pdb file.

I tried those coordinates and managed to create a cloud of 3d points but it is not anywhere near the structure that is found on pdb site.

So wondering how the structure is created on the pdb website?

# Answer
> 4 votes
JSmol sounds like what you want - it'll let you display 3D structures from PDB files onto webpages using JavaScript. If you just want to view the structure rather than display them on webpages, there's many programs that exist already like Swiss-PdbViewer and Pymol. There's even web-interfaces that let you upload a PDB file and they display it for you like Jena3D Viewer.

# Answer
> 4 votes
If you want to share a PDB file, you could try Michelanglo. It allows you to upload a PDB file (among other things) edit a description panel (which can feature special links that control the protein view and representation) and share the link with whomever —without needing any installations on either side. Here is the page of shared protein: gallery (disclaimer: I am the admin).

If you are interested in implementing a molecular viewer on your site, I'd suggest NGL. The PDB website uses NGL too, but is recently switching to MOL\* ("molstar"), which is a new development, so less documented, but fixes some corner cases and is marginally faster and is a collaboration with a PDBe. NGL does not depend on any technology that will be depracated mind you, so it actually is a very solid choice for a JS molecular viewer.

The PDB format is a working format, not a deposition format, as is mmCIF. As a result NGL and Mol\* have some issues when it comes to PDB files, such as lack of secondary structure records (HELIX/SHEET) results in anaemic sheets, fancier linkages (LINK/SSBOND reconds) are not handled only CONECT records.

# Answer
> 2 votes
As others said, there already exist javascript libraries to visualize a protein. The most complete and recent one would probably be Mol\*, which is what is used on RCSB PDB website. Mol\* actually works using NGL, a WebGL engine. But probably the oldest library used on RCSB PDB to visualize proteins would be JSMol. You can actually switch between these viewers on RCSB PDB 3D view using "Select a different viewer" in the bottom right (try it out on 1UBQ)

Now to address the first exercise: there are several things going on after getting the atom coordinates. A PDB file may or may not contain a CONECT part, where bonds between atoms will be clearly stated (contrarily to a CIF file which I believe always contain the bonds). Henceforth these packages may or may not have to guess atom bonds. Given that a PDB file also contains the residues of each atom it's fairly easy to do, less so when dealing with small molecules.

However what you see most of the time is a cartoon representation of proteins, where each secondary structure is represented differently. The library will compute the secondary structure of each residue based on the dihedral angles of the backbone, then represent them. As of why they look like this, check the work of Irving Geis and Jane Richardson.

So if you want to have a website showing protein structures I would strongly advise to use an already existing library, although the exercise of coding your own can always be interesting. If you want to visualize structures on your own computer without using a web interface I can only recommend PyMOL, VMD and ChimeraX.

# Answer
> 2 votes
I am not good as you are and needed to use a parser for the pdb file

see https://mmcif.wwpdb.org/docs/software-resources.html:

```
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  7 17:40:03 2020

@author: Pietro

"""

import matplotlib.pyplot as plt  #matplotlib

from pdbx.reader.PdbxReader import PdbxReader #PDB file parser

import numpy as np #numpy

openz = open('./1msc.cif')

print('opened')

# Initialize a PdbxReader object with the input file handle

pRd = PdbxReader(openz)

# Initialize a list to be propagated with DataContainer (and/or DefinitionContainer) objects (of the DataContainer class, which inherits from ContainerBase) parsed from the CIF file, where data blocks map to DataContainer objects

data = []

print(data)

# Call the read(self, containerList) method with your list

pRd.read(data)

print ('data ###################################################################\n' , data)

print('data lenght' , len(data))

# Your list is now propagated with one or more DataContainer objects, which represent data blocks. To get the first data block, just use list notation:

block = data[0]

print(' block ##############################################################\n',block)

pippo=block.getObjNameList()

print(pippo)

# To retrieve a category object, use the getObj(self, name) method

atomsite = block.getObj("atom_site")

print('atom site : ', atomsite)

atom = atomsite.getValue("group_PDB",0)

# print('atom sites id :', atom)

atomlist=[]

i=0
while True:
    atom =  atomsite.getValue("group_PDB",i)
    atomid = atomsite.getValue('label_atom_id',i)
    if atom =='ATOM':
        # print(atomsite.getValue("group_PDB",i))
        if atomid == 'CA':                         #select only carbon alpha (CA) to work with

            atomlist.append([float(atomsite.getValue('Cartn_x',i)),\
                           float(atomsite.getValue('Cartn_y',i)),float(atomsite.getValue('Cartn_z',i))])

        i +=1
    else:
        print('loop ends')
        False
        break

atomlist=np.array(atomlist)

print('end')

print(atomlist)

print('lunghezza : ', len(atomlist))

fig = plt.figure()
ax = plt.axes(projection="3d")

for p in atomlist:       #print CA as dots
    ax.scatter3D(p[0], p[1], p[2])
    print(p)

for i in range(0,len(atomlist)-1):   #connect with a line CA n with n+1
    x,y,z = [atomlist[i+1][0],atomlist[i][0]],[atomlist[i+1][1],atomlist[i][1]],[atomlist[i+1][2],atomlist[i][2]]
    ax.plot(x,y,z, color='red')

plt.show()

```

my input is https://files.rcsb.org/download/1MSC.cif

Result is just the CAs (carbon alpha) shown as dots and connected to form the main chain:

Didn't manage how to show the dimer though

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


# Question
Title: Number of sing-copy genes identified by orthofinder is far less than BUSCO assessment?
Dear Orthofinder users and BUSCO users,

Recently, I have assessed 30 genome annotataions `using BUSCO with lineage` dataset `ascomycota_odb10`, I got 1100+ complete and single-copy genes. But the `orthofinder` only identified 300+ single-copy genes. Aneone else have encountered this similar problem?

Thanks,

Yanpeng

# Answer
> 1 votes
It's hard to say for sure without seeing the exact commands you used. However, Orthofinder identifies/defines single-copy orthologues as those genes present in a single copy in all species being compared in the analysis. If a BUSCO orthologue is present in more than one copy in one of the species being compared (could be multi-copy, could be fragmented etc.), then it won't be recovered as single copy by Orthofinder.

---
Tags: orthologues, busco, orthofinder
---


title: What is the meaning of split read?

text:
# Question
Title: What is the meaning of split read?
I want to use rna seq data to later perform functional tests on fusion genes. so before that I need to filter the "best results" (of rnaseq) for deciding which candidates I actually want to try in the lab. one criterium is the split read. (among coverage, discordant mate and confidence)

In some cases the split read (or discordant mate/coverage) of one of the two genes that fusioned is "0".

Does that automatically mean that the gene is no longer a candidate for functional testing?

Unfortunately I only got the results in an excel data sheet and don't know which tool was used to gain data.

# Answer
For discordant mapped reads, a position of 0 for the pair of a mapped read usually means that the pair couldn't be mapped anywhere. This would be expected for fusion events because some reads could originate from across the fusion site (i.e. not being fully contained in either "unfused" gene), leading to a sequence that doesn't match anything in the reference.
> 1 votes

---
Tags: rna-seq, sequence-alignment, read-mapping, rna-alignment
---


# Question
Title: Timeout when downloading the ncbi nr blast database
I am experiencing timeout problems when downloading the NCBI nr preformatted blast database using the update\_blastdb script (version 504861).

I run the script with the following paramters

```
update_blastdb --decompress --passive --verbose nr

```

and I get the following error message (in verbose mode)

```
Downloading nr (45 volumes) ...
Downloading nr.00.tar.gz...Net::FTP=GLOB(0x5610fb59b8f8)>>> PASV
Net::FTP=GLOB(0x5610fb59b8f8)<<< 227 Entering Passive Mode (165,112,9,229,195,144).
Net::FTP=GLOB(0x5610fb59b8f8)>>> RETR nr.00.tar.gz
Net::FTP=GLOB(0x5610fb59b8f8)<<< 150 Opening BINARY mode data connection for nr.00.tar.gz (18745730730 bytes)
Net::FTP: Net::Cmd::getline(): timeout at /usr/share/perl/5.26/Net/FTP/dataconn.pm line 82.
Unable to close datastream at /usr/bin/update_blastdb line 202.
Net::FTP=GLOB(0x5610fb59b8f8)>>> PASV
Net::FTP: Net::Cmd::getline(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 203.
Failed to download nr.00.tar.gz.md5!
Net::FTP: Net::Cmd::_is_closed(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 101.
Net::FTP: Net::Cmd::_is_closed(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 101.

```

The timeout happens after ~35 minutes and a file that is approximately 18GB big is being downloaded, which matches the expected filesize. The checksum file (nr.00.tar.gz.md5) is not downloaded. So I'm not sure which of the two files is actually the problem.

I tested downloading the nt database and everything seems to work fine, so I don't think the script is the problem. For comparison, this is the nt download output for the first file

```
Downloading nt.00.tar.gz...Net::FTP=GLOB(0x562575a73168)>>> PASV
Net::FTP=GLOB(0x562575a73168)<<< 227 Entering Passive Mode (165,112,9,229,196,51).
Net::FTP=GLOB(0x562575a73168)>>> RETR nt.00.tar.gz
Net::FTP=GLOB(0x562575a73168)<<< 150 Opening BINARY mode data connection for nt.00.tar.gz (4065912989 bytes)
Net::FTP=GLOB(0x562575a73168)<<< 226 Transfer complete
Net::FTP=GLOB(0x562575a73168)>>> PASV
Net::FTP=GLOB(0x562575a73168)<<< 227 Entering Passive Mode (165,112,9,229,195,107).
Net::FTP=GLOB(0x562575a73168)>>> RETR nt.00.tar.gz.md5
Net::FTP=GLOB(0x562575a73168)<<< 150 Opening BINARY mode data connection for nt.00.tar.gz.md5 (47 bytes)
Net::FTP=GLOB(0x562575a73168)<<< 226 Transfer complete

```

Any help would be appreciated.

# Answer
> 1 votes
The answer of the support stated that it is likely that the size of the first download is the problem and that I should try using another tool like rsync. Because I'm unfamiliar with rsync I decided to write a python script that did the job.

```
import urllib, tarfile, json

base_url = "https://ftp.ncbi.nlm.nih.gov/blast/db"

# download manifest file to get the filenames
manifest_file = 'blastdb-manifest.json'
urllib.request.urlretrieve(f"{base_url}/{manifest_file}", manifest_file)

with open(manifest_file) as f:
    manifest_data = json.load(f)

# download everything
for file in manifest_data['nr']['files']:
    # download checksum file
    checksum_file = f"{file}.md5"
    urllib.request.urlretrieve(f"{base_url}/{checksum_file}", checksum_file)

    # download the archive file
    urllib.request.urlretrieve(f"{base_url}/{file}", file)

    # check that the checksums match
    calculated_checksum = get_md5_for_file(file)

    with open(checksum_file) as f:
        for line in f:
            line = line.strip()
            if line.endswith(file):
                checksum = re.split("\s+", line)[0]
                if checksum == calculated_checksum:
                    checksum_matches = True
                else:
                    raise Exception(f"Checksum doesn't match expected for file {file}")

    # unpack the archive
    tar = tarfile.open(archive_file)
    tar.extractall(file)
    tar.close()

```

The function to calculate the md5 checksum is the following (borrowed from this question)

```
def get_md5_for_file(file):
    hash_md5 = hashlib.md5()
    with open(file, 'rb') as f:
        for chunk in iter(lambda: f.read(4096), b""):
            hash_md5.update(chunk)

    return hash_md5.hexdigest()

```

Note that this is not the exact script I use, so I can't guarantee it works exactly like this, but I tried to include all the important stuff.

In my version I include checks to whether my existing blast database needs to be updated (I keep the manifest files and compare the manifest\['nr'\]\['last\_updated'\] values) and I also try to skip already unpacked archives while iterating (check if the nr.xx.phd file exists) to save time in case the script fails due to problems in the download.

# Answer
> 1 votes
Same here. I used axel to download it a week ago, but fail to get finished. Also, the downloaded taxa zipped file (60Mb) doesn't match md5sum. Apparently, they update the FTP files every 30 minutes or 1 hour.

---
Tags: blast, phylogenetics
---


# Question
Title: How to detect bedtools version used by pybedtools (in order to correctly preserve strand information when merging gtf records)
I have issues working with pybedtools due to strand information being lost after extracting and merging transcript coordinates from a gtf file.

It seems that the solution to preserve strand information depends on the version of bedtools, so I would like to have my code detect bedtools version and adapt its behaviour accordingly.

On a machine with bedtools v2.27.0:

```
ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).saveas()
print(my_genes)

```

```
I       ensembl transcript      4119    10230   .       -       .       gene_biotype "protein_coding";gene_id "WBGene00022277";gene_name "homt-1";gene_source "ensembl";gene_version "1";p_id "P18791";transcript_biotype "protein_coding";transcript_id "Y74C9A.3";transcript_name "Y74C9A.3";transcript_source "ensembl";transcript_version "1";tss_id "TSS21161";
X       ensembl transcript      17714968        17718714        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8395";transcript_biotype "protein_coding";transcript_id "cTel55X.1b";transcript_name "cTel55X.1b";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";
X       ensembl transcript      17714968        17718717        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8874";transcript_biotype "protein_coding";transcript_id "cTel55X.1a";transcript_name "cTel55X.1a";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";

```

```
merged = my_genes.sort().merge(s=True).saveas()
print(merged)

```

```
I       4118    10230
X       17714967        17718717

```

(**Missing strand information**)

```
merged = my_genes.sort().merge(s=True, c=7, o="distinct").saveas()
print(merged)

```

```
I       4118    10230   -
X       17714967        17718717        +

```

On a machine with bedtools v2.26.0:

```
ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).saveas()
print(my_genes)

```

```
I       ensembl transcript      4119    10230   .       -       .       gene_biotype "protein_coding";gene_id "WBGene00022277";gene_name "homt-1";gene_source "ensembl";gene_version "1";p_id "P18791";transcript_biotype "protein_coding";transcript_id "Y74C9A.3";transcript_name "Y74C9A.3";transcript_
source "ensembl";transcript_version "1";tss_id "TSS21161";
X       ensembl transcript      17714968        17718714        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8395";transcript_biotype "protein_coding";transcript_id "cTel55X.1b";transcript_name "c
Tel55X.1b";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";
X       ensembl transcript      17714968        17718717        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8874";transcript_biotype "protein_coding";transcript_id "cTel55X.1a";transcript_name "c
Tel55X.1a";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";

```

```
merged = my_genes.sort().merge(s=True).saveas()
print(merged)

```

```
I       4118    10230   -
X       17714967        17718717        +

```

```
merged = my_genes.sort().merge(s=True, c=7, o="distinct").saveas()
print(merged)

```

```
I       4118    10230   -       -
X       17714967        17718717        +       +

```

(**Duplicate strand information**)

Given that I work on linux systems, **how can I detect bedtools version from within python**?

**Is there a better way** to do what I want?

(Ideally, I would like to have my resulting bed consist in chrom, start, stop, name, score, strand fields, where the score would just be "." and the name would be the content of `itvl["gene_id"]`)

---

Edit (05/05/2021): Regarding my last comment above, I found out that it is easy to convert gff into bed using pybedtools.featurefuncs.gff2bed (tested with bedtools 2.30):

```
from pybedtools import BedTool
from pybedtools.featurefuncs import gff2bed

ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).each(gff2bed, "gene_id").saveas()
print(my_genes)

```

```
I       4118    10230   WBGene00022277  .       -
X       17714967        17718714        WBGene00007068  .       +
X       17714967        17718717        WBGene00007068  .       +

```

```
merged = my_genes.sort().merge(s=True, c=[4, 5, 6], o="distinct").saveas()
print(merged)

```

```
I       4118    10230   WBGene00022277  .       -
X       17714967        17718717        WBGene00007068  .       +

```

This ensures proper bed format.

However, with bedtools 2.26, `s=True` makes strand information appear as a 4th column, so it is more complicated to obtain correct column order after merging.

# Answer
> 3 votes
pybedtools assumes that bedtools is in your path and bedtools itself will return the version with `bedtools --version`. So:

```
import subprocess
subprocess.check_output(['bedtools', '--version'], text=True).strip().split()[1]

```

For me that returns `'v2.30.0'`.

---
Tags: python, bedtools, pybedtools, strandedness
---


# Question
Title: Use Htslib to extract specific region of vcf file
I am working in c++ using the htslib library to parse different vcf files as well as fasta files.

My goal is to read in a specific region from a fastafile and then extract all the VCF sites within that region from a seperate vcf file, in this case only containing chr14. I can extract the fasta region in several ways, one of which is to use the function faidx\_fetch\_seq() and define a region of chr14:20000000-22000000.

Then for the bcf files i am able to iterate through the entire file and extracting the different information regarding the reference, alternative allele and so on. With a simple example shown below.

```
void faVCF(const char* fastafile,FILE *fp){
  htsFile *test_bcf = NULL;
  bcf_hdr_t *test_header = NULL;

  test_bcf = bcf_open("chr14.vcf", "r");
  test_header = bcf_hdr_read(test_bcf);
  bcf_hdr_set_samples(test_header, "20000000-22000000", 0);
  bcf1_t *test_record = bcf_init();
  while(bcf_read(test_bcf, test_header, test_record) == 0) {
              std::cout << bcf_hdr_id2name(test_header, test_record->rid) << "\t" <<
                          test_record->pos << "\t" <<
                          test_record->n_allele << "\t" <<
                          std::endl;
          }
}

```

However, i doesn't work to extract all the variant sites in a specific region. I tried bcf\_hdr\_set\_samples with the argument "20000000-22000000", to extract that given region from my chromosome. As far as I understand this should be able to subsample both samples and region queries using the header information, and then upon initiation extract that region, however this doesnt work.

So my question is: Can anyone help me with what function I can use to extract all the vcf sites from the region chr14:20000000-22000000.?

# Answer
For iterating over reads in a region in a single BCF file, you want to use the iterator interface: `bcf_itr_queryi()` or `bcf_itr_querys()`, together with `bcf_itr_next()`. For a VCF file, it appears you need to use the similar TBX iterator interface: `tbx_itr_queryi()` or `tbx_itr_querys()`, together with `tbx_itr_next()`. See HTSlib's *tabix.c* for reasonably straightforward examples of both.

Alternatively you could use the routines in *htslib/synced\_bcf\_reader.h*, which encapsulates both and is what bcftools mostly uses.
> 4 votes

---
Tags: c++, htslib
---


# Question
Title: Finding 2000 bp sequences that occur more than once in the genome
I am looking for a software that can identify long sequences (1500-2000 bp) that occur more than once in my genome of interest. A k-mer counter like KMC could have worked, but KMC has a max limit of 256 bp. If I can club contiguous 256-mers from KMC, that would work too, but I don't know how to do that. Any ideas?

# Answer
> 3 votes
MUMmer includes a `repeat-match` tool which can do exactly that:

```
$ repeat-match
USAGE:  repeat-match  [options]  <genome-file>

Find all maximal exact matches in <genome-file>

Options:
 -E    Use exhaustive (slow) search to find matches
 -f    Forward strand only, don't use reverse complement
 -n #  Set minimum exact match length to #
 -t    Only output tandem repeats
 -V #  Set level of verbose (debugging) printing to #

```

Note that a "maximal exact match" is just a maximal repeat.

---
Tags: sequence-alignment, k-mer, repeat-elements, sequence-homology
---


# Question
Title: DESeq2 multiple treatments, multiple time points, multiple cell lines
*Note: this question has also been asked on Bioconductor Support*

I know there are a lot of questions asking similar things here on this forum and I checked the vingette and did a lot of other research but I still have a hard time wrapping my head around it.

In my position, data is unfortunately just put in front of me and I'm asked for pretty pictures. Even making it clear to the wetlab people in the lab that we need replicates is an issue.

I have multiple cell lines, multiple time points and multiple treatments:

Cell lines: CL1, CL2, CL3 Time points: 6h, 24h Treatments: T1, T2, T3, Control

For each cell line and each time point, there are 3 different treatments plus a control. 3 replicates for each sample -\> 72 samples

What I want is actually quite simple; I want to measure control vs each treatment at each time point in each cell line. No testing across time points or across cell lines.

My first thought was to separate the data into 6 different data sets (CL1, 6h | CL1, 24h | CL2, 6h | CL2 24h | CL3, 6h | CL3, 24h) and simply do the DE analysis separately, but I read that this is not the way to go. In the end it would also be nice to get a normalized count matrix with all data normalized together, for PCAs and similar.

I hope I made it understandable. How do I design the DESeqDataSet?

# Answer
As Mike Love alluded to on the Bioconductor pages, it's best to have specialised help for datasets that are this complex; doing analysis on anything which has multiple compounding factors is not simple. It's also a good idea to talk through the analysis *before* doing sequencing, so you can discuss things like having six replicates per sample (including controls).

In any case, what I did for something similar to this (eventually, after lots of discussion with biologists about possible models) was to combine all the non-tested factor columns into a single factor. Something like this:

```
meta.combined.df <- meta.df %>%
  mutate(lineTime = paste0(cellLine, "_", timePoint))

```

And then subset the data based on this statistic:

```
counts.sub.mat <-
   counts.sub.mat[,meta.df$lineTime == desired.lineTime]
meta.sub.df <-
   subset(meta.combined.df, lineTime == desired.lineTime)

```

Then use `~ Treatment` as the design formula (assuming "Treatment" is a field containing either the treatment type or "control").

But... this is explicitly treating each line/time point as a separate experiment, which may not be what you want. There are complications associated with batch effects and outliers that need to be talked through, because you could end up identifying "differential expression" that is not reproducible.
> 1 votes

# Answer
To compare a subset of samples to another subset of samples, make one column of column data which will distinguish them all the way you want (so probably join cell line, treatment, timepoint), and make that column your design, and then when you call results, you specify the pair you want to contrast, one pair at a time.
> 1 votes

---
Tags: deseq2
---


# Question
Title: Problems about value log2(IP/Input) less than zero in ChIP-seq?
I have a question about the normalization for ChIP-seq. I used CPM to normalize my bam files of each IP and Input. Then I calculate the coverage of gene bodies for all genes on the genome. I have WT and different mutants. Take H3K4me1 as an example, I first calculate the matrix over gene bodies for H3K4me1 IP and Input. Then I used the matrix of H3K4me1 divided by the matrix of input. then log2(H3K4me1/Input) and plot this value over entire gene bodies. I found all the regions are lower than 0 and it seems that my IP does not work. In fact, this is not the case since I have checked the antibody's efficiency via ChIP-qPCR and it works well. Moreover, the distribution of H3K4me1 is consistent with the published data. But the problem is the profile of log2(H3K4me1/Input) is lower than 0. I don't know how to deal with such problems? I would appreciate it if someone can give me some suggestions on that. Thanks!

# Answer
Seconding what Devon says:

1. Use a proper normalization technique such as TMM from `edgeR`. For applying this to bigwig files you could follow my tutorial at biostars.org.

My two cents beyond that:

2. I would not use enrichment over input as input is hard to normalize and can be sparse depending on sequencing depth. Therefore, you might get biased enrichments that do not well reflect reality, especially if counts in IP or input (of both) are low. I personally only use input for peak calling, and for plots like this then compare IPs among each other.

So the strategy would be to first normalize your bigwigs (I guess you use these for the plots?) using the strategy as linked above, and then use these robust CPMs to make the plots without any input samples, only using the IPs. You can use MA-plots to explore whether samples are properly normalized. I made a little code suggestion here in this biostars post on how to do this starting basically from a matrix of raw counts using `edgeR`. This is basically (from the TMM normalization standpoint) the same as in the first link, but then doing pairwise MA-plots to see whether TMM properly centers the majority of points along y-axis being somewhat zero. Normalization is key in any experiment, therefore I always recommend to use these kinds of plots for diagnostics.
> 1 votes

# Answer
This is a common issue caused by the normalization method not being particularly robust (most likely it's using the mean signal for normalization, which will be a bit upward skewed when you have regions of high enrichment in the ChIP). The basic ways to fix this are:

* Use a more robust normalization (e.g., the median signal, or even just manually select un-enriched regions to calculate the normalization factors)
* Manually tweak the normalization factors so the background (i.e., log2 ratio in regions where you don't expect any enrichment) is 1. That will generally make the regions you're plotting more biologically realistic.

Note that there can be a negative enrichment in some areas, due to things like exclusion of histones, so don't be overly zealous if you manually tweak the normalization factors.
> 0 votes

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


# Question
Title: Calling for one specific SNP from multiple sequencing runs
Just learning the basics of bioinformatics and bash shell scripting.

Is there a correct way to search for a specific allele with one or even multiple sequencing runs? I am searching for prevalence of the C2416T allele in the SARS-CoV-2 genome within the ORF1a RNA polymerase gene as a mini project.

I downloaded a bunch of SARS-COV-2 sequencing runs, I think I need to convert them to bam file type than I am not sure how to call for one specific SNP, especially in a loop for 5-7 runs. Can this be done with bash?

I was looking into variant calling from a bam file type using vcf and I can't find a way to call for one specific snp on multiple bam files at once.

Ideally I would want to filter the variant out of the .bcf files or use grep since I sort of understand those the most and I’m not working with a large amount of data or a human genome.

# Answer
> 4 votes
Getting SNPs from whole-genome data is usually done via mapping sequencing reads on the reference (bwa-mem or bowtie are two popular mappers) and then using one of variant calling tools (like freebayes, samtools, or GATK). However, I would also look up if there is something virus-specific.

According to this answer, you could get variant calls of a specific region using samtools if you give it a bed file that contains the region of interest (three columns `chr start stop` format). For example using human reference

```
samtools mpileup -uDl regions.bed -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - > RAL_samtools.raw.bcf

```

You could also specify a bunch of files, not just two and it would do the variant calling on them all.

# Answer
> 2 votes
*Answer from @atpoint \[transferred from comment\]:*

Alternatively, you can use `samtools view` to only retain reads that overlap the region you are interested in and then feed these into the mpileup. IIRC there was something with that `-l` option in mpileup that caused trouble but I do not remember what it was. I think it still called all variants based on all reads and then simply filtered for the target regions in `-l`. The `view` approach would probably be much faster. Just keep that in mind as an alternative.

---
Tags: variant-calling, bash, covid-19, sars-cov-2
---


# Question
Title: Base abbreviations other than ACTGU in sequence file
**Short question:**
where can I find what `F` and `J` and `<` and `-` and `7` mean in a Sequence file?

**Longer version:**
I was looking at this sequence file I got from NCBI (it's an sra file? if that is the correct terminology) and I used the `sratoolkit` to print out a few `spot`s

looking at them I expected to see ACGU (this was meant to be an RNA sequence), but what I see are like

> AAFFFJJJFJJJJJJJJJJJJJJJJJJJJJFJJF\<JJ77AJJJJ-JJJJJFF7FA--F7FAJJ7\<AJ-FJA\<J\<JJFA\<\<77F-JJJJJJF\<JF\<7-FF-A-----A-77JJ7FJJF-FFFAJ\<\<-FJJFAJAFFAFFFFF\<7A7JFJJ7AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJFJJJJJJJJJJJJAA\<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJAAFFJJFJJJJJJJJJJJJJJFJ-FJF

Where can I find a reference to help me "understand" this?
I mean I assume `-` means missed or unknown? or is it meant to be already aligned and these are deletions?

I appreciate any pointer to get started. as you must have observed I am a total noob.

I got the file from: https://www.ncbi.nlm.nih.gov/sra/SRX8156420\[accn\]

Then I did a

```
 ./fastq-dump -X 5 -Z ERR4557575

```

Follow-up question:

* Where do you look for the quality scores ranges and such, I guess as @winnie2k suggested I should convert this to FASTQ and make my life easier, but still wanted to understand the file better.

# Answer
> 3 votes
I usually convert SRA files to FASTQ before doing any further processing, but the tooling for converting SRA files is cumbersome.

To simplify things, I recommend you download the FASTQ file directly. The ENA currently has a nicer interface for doing these things. Here is a link to download the FASTQ files for the SRA sequencing run associated with the SRA experiment you are describing: https://www.ebi.ac.uk/ena/browser/view/SRR11588968

A description of the FASTQ format can be found here: https://en.wikipedia.org/wiki/FASTQ\_format

---
Tags: rna-seq
---


title: Shannon's Entropy and DNA

text:
# Question
Title: Shannon's Entropy and DNA
I have been working with Shannon's entropy (SE) and DNA lately, and referring to the formula and concept of SE, I just wondered whether one should in the case of DNA use each nucleotide as having a probability to show up as 0,25 or if it was better to use the probability based on the frequency in the genome (let's say the human genome).

In the case of the latter being the better, is it interesting to have those frequencies at the scale of the genome, the chromosome, the sub segment ?

I'm interested in calculating the entropy relative to a short window (the length of an Illumina read, about 125-150bases), using the proportion of each base on a given sequence length. The idea is "how much information is given" in this sequence, and by information I mean how much different bases do I have.

It might that I'm missing some basic mathematics knowledge here, don't hesitate to point it out so I can understand it better

cheers,

# Answer
> 1 votes
I think that genome-wide probabilities (i.e. not 0.25 for each base) should be fine. For a given segment of DNA, the chance that the sequence doesn't appear anywhere else in the genome is fairly high.

There is lots of patchy non-randomness (e.g. at centromeres, STRs, and other tandem repeat regions), but there's not much point in modelling that because it's not random.

---
Tags: phylogenetics
---


title: 10bp sliding window quality score

text:
# Question
Title: 10bp sliding window quality score
I am having trouble creating a 10 basepair sliding window that goes through `seq` and trim them all 10 basepairs that have an average quality score (`qual`) less than 15. I am trying to approach this by reading in the file and seperating the sequences and quality scores, then read in the 10 basepair window in the sequence.

<pre><code>#!/usr/bin/envpython
with open('fastaQ.fq') as f:
    lines = [line.rstrip() for line in f]
    seq= []
    qual= []
    tenFrame= []

    #Seperates the sequences and the quality score of the file
    for i in range(len(lines)):
        if lines[i].startswith('@'):
            seq.append(lines[i+1])
            qual.append(lines[i+3])
            i=i+4

    #
    for i in range(0,len(seq)):
        tenFrame.append(seq[i][i:i+10])
        if len(tenFrame)%10 == 0:
            print(tenFrame)```
</code></pre>

# Answer
I'm not sure what you mean by "trim". The code below will discard any reads that have a 10bp window with a mean phred score below your cutoff. I'll ammend the post to also include a solution that "trims" basepairs on both ends until the window achieves an average score of above the cutoff. It would be ideal if you could provide a minimal reproducible example that shows what your code does and what you'd like it to do.

As far as the code goes, I'd first suggest using a library like Biopython to parse the Fastq file. It will make your life a lot easier. For your code, it looks like `seq` actually contains all the sequences, so `len(seq)` is the number of reads in the file. However, when you iterate over `i in range(len(seq))`, it looks like you're using `i` to both index which read to iterate over as well as get the 10 base pair window. You need an additional index variable to obtain the window:

```
for window_start in range(len(seq[i]) - window_len) + 1)
    window = seq[i][window_start:window_start + window_len]

```

Additionally, some pythonic tips for the future:

* If a range starts from 0, don't specify it as an arugment i.e. use `range(x)` as opposed to `range(0, x)`
* Try to avoid anything along the lines of `range(len(x))`. In this case, you could just iterate over the sequences directly: `for s in seq:` or, if you wanted the index, you could enumerate the sequences via `for seq_idx, s in enumerate(seq)`
* Use good variable names. `seq` implies a single sequence, yet is a list of sequences. Same with `qual`

Here is some code that will retain only sequences that have all windows with a passing mean phred score

```
from Bio import SeqIO

read_file = "/file/here.fq"
w = 10
cutoff = 15
valid_seqs = []
with open(fq_file) as fq_fd:
    seqs = SeqIO.parse(fq_fd, "fastq")
    for seq in seqs:
        scores = seq.letter_annotations["phred_quality"]
        mean_scores = (np.mean(scores[i:i+w]) for i in range(len(scores) - w))
        if any(mean_score < cutoff for mean_score in mean_scores):
            print(f"{seq.id} does not pass the filter!")
            continue
        valid_seqs.append(seq)

```
> 1 votes

---
Tags: python
---


title: Cellranger gives error

text:
# Question
Title: Cellranger gives error
I am trying to run cellranger but I get fastq permission denied error

```
fi1d18@RBGO-Server2:/opt$ cellranger-atac count --id test --reference refdata-cellranger-arc-GRCh38-2020-A-2.0.0 --fastqs /data/fi1d18/Downloads/FASTQ_scATAC/ --chemistry ARC-v1
error: couldn't open MRO file
    Caused by: Permission denied (os error 13)
fi1d18@RBGO-Server2:/opt$

```

I have added the path by

```
export PATH=/opt/cellranger-atac-2.0.0:$PATH

```

And by vim I have also added the path

```
fi1d18@RBGO-Server2:/opt$ vim ~/.bashrc

```

Any help please?

# Answer
Two of those command line parameters should be accessible locations (see here). Could you please add the output of the following commands:

`ls -l /data/fi1d18/Downloads/FASTQ_scATAC/ | head`

and

`ls -l refdata-cellranger-arc-GRCh38-2020-A-2.0.0 | head`

I suspect that the refdata link should be an absolute path, rather than a relative path, i.e. something like this:

```
cellranger-atac count --id test --reference /opt/refdata-cellranger-arc-GRCh38-2020-A-2.0.0 --fastqs /data/fi1d18/Downloads/FASTQ_scATAC/ --chemistry ARC-v1

```
> 1 votes

---
Tags: single-cell, bash, linux, 10x-genomics, cellranger
---


title: Chimera Alignments

text:
# Question
Title: Chimera Alignments
I have a structure with two subunits. I am trying to show movement of the C-terminal subunit upon ligand binding by superposition with another structure from the same strain in the apo form.

I want to superpose the N-terminal subunits (A), note the root mean square of the atomic deviations (RMSD), then apply the resulting matrix to the superposition of both subunits of the molecule with both subunits of the other (B), then I want to superpose on the C-terminal subunits (C), and then measure the RMSD between (C) and the C-terminal subunits from (B).

Do you know how can I do this in Chimera? I know how to superpose molecules and calculate the RMSD in Chimera, but I do not know how to apply the resulting matrix from (A) to (B). I tried to save the aligned structures from (A) as a PDB file to check and see if there was a matrix in the file, but an error message came up saying that I needed to save the file as `$name` or `$number`. I tried to save the file with `$` in front of the name, but it didn't work.

How can I do this correctly?

# Answer
> 1 votes
*Answer from @terdon converted from comment:*

The `$` are often used to signify that something is a variable. So I guess they mean you need to use a name that you've already used. I mean, the same name as used for something else in your pipeline, so that then the program will see the file with the same name and know to apply it. Sorry, I don't have any experience with chimera, so I have no idea if what I'm suggesting is reasonable. It's just what I get from the `$` you mentioned.

*Answer from @marcin converted from comment:*

this question was answered on chimera-users:

http://plato.cgl.ucsf.edu/pipermail/chimera-users/2017-October/014045.html

*Answer from @joe-healey converted from comment:*

You are probably looking for the “Morph Conformations” function within Chimera. It’ll move and superimpose 2 structures on to one another.

---
Tags: sequence-alignment, software-recommendation, protein-structure
---


title: Help with 1D^2 library shearing

text:
# Question
Title: Help with 1D^2 library shearing
I've done my 1st trial of 1D^2 whole genome sequencing using LSK-SQK308 kit with R9.5 flowcell. Without doing library fragmentation, I encountered the unexpected library shearing during sequencing.

But I didn't see significant portion of low-molecular-weight fragments in the QC gel photo before and after library preparation.

Anyone using nanopore has encountered this before? Any explanation? And how to avoid it? I could not find evidence of the problem in my QC steps...

Thanks in advance!

# Answer
Short template DNA is a big problem with nanopore sequencing; it doesn't easily show up on a gel if it's well distributed across different lengths (as appears to be the case here from the read length plot), and sucks up both adapter and sequencing capacity.

The usual way to deal with it (for reads in the size ranges seen here) is to reduce the bead concentration when doing magnetic bead cleanup. Try something like 0.6X beads (e.g. 30μl of Ampure XP beads added to 50μl template DNA solution), which should still allow most of the long template sequences to bind, but leave shorter reads in solution.
> 1 votes

---
Tags: nanopore, quality-control, 1d2-reads
---


# Question
Title: Is there a modern alignment tool tailored for transmembrane regions?
I am looking for a project or tool that allows programmatic pairwise alignments of proteins but that takes care with transmembrane regions of proteins. TM regions are traditionally too information poor for many tools to align accurately (Wong et al, 2010).

* TM-aligner doesn't appear to have a program to download, and is only available via webserver. On the otherhand,
* PRALINE has some optimisation for TM regions, however was released in 2008 and since then there has been a boom in TM knowledge from structures so does not benefit from that.Currently I am using PSI-BLAST. However,
* PSI-Blast uses SEG to mask low complexity regions. This will most likely mean that TMHs are masked making the alignment of those regions fairly unreliable.

Primarily I am interesting in sequence alignments. But if anyone knows of a TM structural alignment, that could be good to know.

I think generally that counting idea works for an ideal helix, but a vast majority of TMHs break convention. I wouldn't expect to see a pattern based aligner that was that general. I'd be interested if anyone knows of a tool that juggles that though!

If I understand correctly, Matteo was saying that TMHs follow a simple patter (HHP/H,HHP/H...~7 times). This is not the case. There are too many examples of TMHs that do not match this pattern to assume this pattern would accurately align TMHs. As @M\__ suggests, the number of TMP structures has not caught up to their soluble counterparts in the database. I've updated the question to explain why PSI-Blast might not be the best option.

# Answer
*Answer from @matteo-ferla converted from comment:*

I assume you mean sequence (MSA), not structure alignment. Despite being a structural bioinformatician, I would slap a big caveat against structure based aligners. TCoffee:Espresso is probably the best structure-based aligner (browser and download and API), but the results are not as solid as some non-structural approaches such as training Psi-blast and aligning with that matrix. I am not familiar with MSA of TM regions, however, so this group may be an exception.

Given that helices are a repeat of 2 hydrophobic residues and 1.5 either polar (core-wards) or also hydrophobic (depending on how the helices pack), a pattern based aligner ought to do well from an evolutionary perspective. If you really want from a structure perspective and don't want to know about the cases where a loss of a turn at the top is compensated by the bottom, say, then TCoffee is definitely worth a try

*Answer from @m converted from comment:*

PSI-blast is a useful tool for deep alignment patterns, albeit it has its critics. I agree the vast majority of in vivo structures lack TM proteins, for different reasons
> 1 votes

---
Tags: sequence-alignment, proteins, covid-19, membrane-protein
---


# Question
Title: Find covariation of different DNA-binding protein binding between two conditions
I have ChIP-seq data of RNAPII and two other factors which we think follow RNAPII during transcription in two different conditions, and I'd like to show that the genes that lose RNAPII also lose the other two factors, and same for gains.

Do you know any methods that can calculate some sort of covariation coefficients?

I know there are methods that calculate covariation of genes between two conditions, but none for different measures on the same gene.

# Answer
*Answer from @stupidwolf converted from comment:*

I would try something like this, for every gene you have a logFC of RNAPII in control versus treatment. Do this for other 2 proteins, and you correlate the logFC. In very rough terms, correlation is the scaled form of covariance, so this will work.

If you have 1000 genes, do correlation between 1000 logFC from RNAPII against 1000 logFC from other proteins.

*Answer from @m converted from comment:*

You use a covariation stat and then perform Pearson's correlation to standardise the response between -1 to 1. An answer of \>0.7 is considered robust and \>0.8 highly robust. ITs an interesting question and you are right covariation is the way to answerr it
> 1 votes

---
Tags: statistics, chip-seq, transcriptome-regulation
---


# Question
Title: How to load output from the kegg API in biopython into a pandas dataframe?
Biopython provides a (bit unintuitive) API to access to the Kyoto Encyclopedia of Genes and Genomes (KEGG). I am trying to make use of it, but the output is quite unhandy as a string. What is the best way to parse the data into a `pandas.DataFrame()`?

```
from Bio.KEGG import REST as Kegg
import io

pd.read_table(kegg.kegg_get('K02545'))

```

# Answer
> 1 votes
*solution from @sören:*

I came up with this solution for my problem:

```
from Bio.KEGG import REST as kegg

def _get_kegg(kegg_id):
    kegg_output = kegg.kegg_get(kegg_id).read()
    results = {}
    for line in kegg_output.split('\n'):
        splits = line.split()
        if not line.startswith(' '):
            if len(splits) > 0:
                key = splits[0]
                value = ' '.join(splits[1:])
                results[key] = value
        else:
            results[key] += ' '.join(splits)
    return pd.DataFrame(results, index=[kegg_id])

_get_kegg_v = np.vectorize(_get_kegg)

def get_kegg_info(kegg_ids):
    if isinstance(kegg_ids, str):
        kegg_ids = [kegg_ids]
    return pd.concat(_get_kegg_v(kegg_ids), sort=False)

```

---
Tags: python, kegg, pandas
---


title: Where did the BIND website go?

text:
# Question
Title: Where did the BIND website go?
I was trying to reach the BIND (**B**iomolecular **I**nteraction **N**etwork **D**atabase) but the website is down.

Where can I find something similar? Has this platform moved somewhere or what?

# Answer
> 0 votes
*Answer from @heathobrien converted from comment:*

looks like that link has been dead for years:

https://www.biostars.org/p/16763/

*Answer from @jason-aller converted from comment:*

The biostars article mentions this:

http://baderlab.org/BINDTranslation

---
Tags: database, interactions
---


# Question
Title: How to import data from cell ranger to R (Seurat)?
I will have some scRNA-seq data. The goal of the experiment will be to see if there is any difference in gene expression between treatment groups using the package Seurat from R.

I have read a tutorial how to do the analyze, but this tutorial does not explain how to import data. I am a total beginner with R, I searched the internet,but the majority of the explanations related to csv/excel/txt... common files. Data are from Cell ranger and spread in 3 files with following file extensions : .tsv and .mtx (barcodes.tsv, genes.tsv and matrix.mtx).

Collaborators ran Cell Ranger and gave these cell ranger output files : `barcodes.tsv`, `genes.tsv` and `matrix.mtx`.

Can someone give me the code to import these kind of data to R ?

Thank you in advance.

update: By "directory" you means "folder" on my computer for example ? So I need to save the 3 files under one folder and then to put name of the folder into brackets after the Read10x function ? I am a total beginner with R (I have used to R commander for statistical tests but that's it...).

update 2: I successed to import data into Seurat with the explanation of StupidWolf. I do not know C language ( I do not know almost any language ..). My system crashed a little further in the analyze after importing, not enough RAM... so now I am learning Linux commands to switch on the cluster.

# Answer
*Answer from StupidWolf, converted from comment:*

There is a function call Read10X in seurat, you can see how to use it with the example in https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/rkit

You place the three files into a directory / folder, then you specify Read10X()... hey how big is your data? you might need to run this on a machine with enough RAM.

You can try downloading the data following instructions here:

https://satijalab.org/seurat/v3.1/pbmc3k\_tutorial.html
> 1 votes

# Answer
Have a look here and the entire tutorial:

https://github.com/hbctraining/scRNA-seq/blob/master/lessons/03\_SC\_quality\_control-setup.md
> 0 votes

# Answer
To convert data from one format into another, you may find some software tool sometimes. But mostly no such tool can be found to satisfy all the requirements that we want.

In such situations, I recommend that you write a C program that reads data from the input file and writes it into the output file. In this way, you have complete flexibility so that you can determine the exact output format and what part of data to ignore or include. To do that, you should be familiar with some concepts in C language such as arrays, loops, strings, and file operations.

For an example of such conversions, you can get ideas from the following page: https://stackoverflow.com/questions/55347322/c-program-to-convert-a-text-file-into-a-csv-file

Of course, such a converter program can be written in other languages too such as C# and Python.

The program can be developed in R language too. You need to be familiar with import functions, loops, strings, and arrays in this language.

We cannot expect that custom tools exist for all our file convertions. Every bioinformatician must be skilled in a standard programming language such as C, C#, Java, Python, or R. Then, he/she will be able to develop any custom program that reads an input file and generates a custom output file.
> -2 votes

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


# Question
Title: Parsing gene expression value for a list of genes from GEO
I am looking for the gene expression data in GEO GSE15543 with the following keywords,

```
(Pancreatic Beta Cells) AND "Homo sapiens"[porgn:__txid9606]

```

I would like to parse the gene expression value of a list of genes from all the samples reported in any particular study.I understand there are two different formats for downloading datasets, SOFT and MINiML.

I have downloaded the soft.gz file and I am using BioConductor to read the contents from soft file using the instructions given here . The file has the content that I 'm looking for , I have checked the contents on the geo2r webpage. I'm still facing challenge in finding out how to fetch out the gene expression values for a list of genes. Any suggestions?

From the files that I download, I wish to parse the gene expression value and the Sample Title and Sample Description so that I can classify the values based on disease state and control state. Any suggestions on which format can be downloaded for the above-mentioned task?

# Answer
GEO provides a GUI for this (just click on the "Analyze with GEO2 button"), you can even check/download the R/Bioconductor code used for the analyses.

Basically you will create groups of samples to compare/contrast and R/Bioconductor will run behind the scenes. Here is a nice tutorial.
> 1 votes

---
Tags: gene-expression, geoquery
---


# Question
Title: Alignment to verify plasmid from forward and reverse Sanger sequencing reads
I am analyzing Sanger sequencing data to verify whether or not the correct protein coding sequence was inserted during cloning.

For each sample, I have both forward and reverse Sanger sequencing reads that cover the coding sequence.

Would it be better if I first align the forward and reverse reads to get a consensus, and then align the consensus to the reference? Or should I just align each read directly to the reference?

# Answer
> 2 votes
Align all reads from a sample at the same time IMO. This will give you all of the information at once, and it should not be difficult. It should be very obvious whenever there's a disagreement between reads or between reads and reference.

I would also recommend inspecting chromatograms directly at all sites of interest (such disagreements, insertion sites, possible errors or mutations).

---
Tags: sequence-alignment
---


title: Imputation of dog genotype

text:
# Question
Title: Imputation of dog genotype
I am an undergrad working on a set of genotype matrix of dogs, it's in csv file format with information on only the each sample as rows, chromosome and its position as columns, and the rest of the matrix is filled with 0,1,2 and missing values.

So I am wondering what would be the best method to impute the genotype with all the data I get? I looked up Impute2, but it seems to require the reference data(recombination map) which I don't have...

Really appreciate whoever answers my question!!

# Answer
> 2 votes
Impute2 is a great, if dated choice. Another great choice is beagle. The advantage of beagle in this case is that it is not as strict about requiring a recombination map. Without a map, results will be poorer, but perhaps good enough.

Do make sure that your sample size including a reference panel is at least 200 in any case.

---
Tags: imputation
---


# Question
Title: Biological datasets for polythetic classification
**Question:** What are some instances of polythetic datasets in biology? In particular, I am looking for a dataset to benchmark a machine learning algorithm optimized via episodic training.

---

**On polythetic classification**

The conventional definition of a conceptual or monothetic class is that its members must possess certain properties in common \[\]. For example, in the MNIST database of handwritten digits, images can be classified into 10 distinct classes (i.e. digits) where all instances of any given class possess the same set of attributes (e.g. zeros are typically written as ellipses).

In contrast, a polythetic class encompasses members that might share a number of characteristics, but none of these attributes is essential for membership of that group or class. For example, one could group green zeros together with red ones in one polythetic class, and red zeros together with green ones in another one, as shown below.

Are there any instances of biological datasets with similar properties? I am particularly thinking of genetics data (e.g. single nucleotide polymorphisms) where perhaps a certain mutation A counteracts the effect of mutation B? Any pointers are welcome.

# Answer
> 1 votes
I may not be understanding exactly what you are talking about, but it seems to me that nearly any set of classes of interest must be polythetic. Most groups of comparable things differ in some characteristics but also have characteristics in common, and none of those characters is 100% identical with membership in the class.

Some examples off the top of my head:

1. Biological species / taxonomic groups (classify according to phenotypic characters)
2. Biological species / taxonomic groups (classify according to molecular sequence e.g. DNA)
3. Cell types (classify according to cell epitopes)
4. Cell types (classify according to transcriptional profile)
5. Ecosystems (classify according to species profile)
6. Disease state (classify according to any one of many molecular profiles)
7. Members of population groups (classify according to genetic variation)

etc.

---
Tags: snp, machine-learning
---


# Question
Title: Is it possible to correct bond order and conectivity problems in an sdf file according to mol2 files?
I was doing a Virtual Screening Search using Pharmer software. The .sdf output file contains a column in mol format and another column with Root Mean Square Deviations. I can see some bond order and connectivity problems in the mol column. I do have a mol2 file that contains the original molecules and I was wondering if I can use these mol2 files, to fix this. Can rdkit work in this case?

# Answer
RDKit has a one line command to fix the bond order (docs)

```
 mol = Chem.MolFromMolFile('wonky.mol')
 template = Chem.MolFromMol2File('template.mol2')
 newMol = AllChem.AssignBondOrdersFromTemplate(template, mol)

```

The problem is that the two need to be identical. I see a carboxylate in your screenshots, so you most likely did the correct thing and docked with the appropriate protonation for pH 7, i.e. as the conjugate base with a negative formal charge. If one does not then RDKit will fail.

Also, RDKit is "not great" with mol2 files. If you have the original SMILES that is better. If the problem is atom names, using a PDB is a better call. If you want to assign atom names, this blog post of mine may help.
> 2 votes

---
Tags: rdkit
---


# Question
Title: Samtools/bcftools calling indels in noisy reads
I have very noisy nanopore reads and am trying to call SNPs/indels.

I'm runinig into some trouble when truing to use the `samtools mpileup | bcftools call` combination.

It seems to incorrectly be calling very long indels where it looks like there is no support. For example:

`samtools mpileup -f ref.fa sample.bam -r Chromosome:198940-198940` produces:

```
Chromosome    198940    C    37    .,.-1A.-1A...-1A.-1A.-1A.,-1a.-1A.,-1a,-1a.,.-1A,,.-1A,,+1a,+1a,.-1A,,-1a.,-1a,-1a,,-1a,-1a.-1A,-1a.    06?5/<:8>/4<?3691465<3354.08:23:11@3C

```

This looks fine... however when I try to call the genotype with bcftools I get a very long indel. `samtools mpileup -ugf ref.fa sample.bam -r Chromosome:198940-198940 -t AD | bcftools call -mv` produces:

```
Chromosome    198940    .    CAAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA    C    225    .    INDEL;IDV=32;IMF=0.5;DP=62;VDB=0.145385;SGB=-0.662043;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,7,2;MQ=60    GT:PL:AD    1/1:255,27,0:0,9

```

Using `htsbox pileup -f ref.fa sample.bam | grep 198940` it shows there is one read with this indels together with many other reads with smaller indels.

```
Chromosome  198940  C   C,C-75AAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA,C-4AAAA,C-2AA,C-1A,T,C+1A,C+1T,C+3GGT 0/4:25,1,1,2,24,5,2,1,1

```

Is there any reason why it would be grouping these reads together into the big indel?

Any idea what is going wrong here?

# Answer
It's likely that the caller is being overly sensitive to deletions, which are the most common error mode for nanopore reads, even three years later in 2021.

What I'd guess is happening is that the caller is seeing a string of consecutive 1bp deletions at high frequency, and treating these as a consensus deletion covering many more bases. Looking at the pileup calls for the surrounding bases may help to diagnose this.

My general advice to people doing variant calling with nanopore is to be suspicious of deletions, especially if they include homopolymers. That usually means trying to confirm any important variants with another technology.

Newer, more accurate, base calling models are being created all the time by ONT, so if you still have the raw FAST5 data, it might be worth it to run through the most recent base caller.
> 1 votes

---
Tags: samtools
---


# Question
Title: Annotatation of transposable elements / transposons / SINE sequences in human
I need to annotate human genome with different types of information, amongst them: transposable elements. I read about repeatmasker, but what I understood is that it estimates the probability of a sequence to be a transposon based on the sequence and not based on curated database. My university does not have access to RepBase (GIRI). I was hoping finding some possibility with biomart but it seems not possible.

What would be a way to access curated annotation of all kinds of transposable elements (SINE, Alu seq..) within a given genome/genome coordinates? My goal is to estimate the link between some measurements we have made on the genome (binned) ( and Shannon entropy is one of them, and probably the less interesting regarding TE). Eventually, I'd like to be able to estimate that measurement X has a correlation with the number of Alu repeats present in the neighbourhood.

What would be the best way to achieve this?

# Answer
> 3 votes
repeatmasker is meant to annotate/mask a genome given a library of known transposable elements. From the project description, I assume you have not done any sequencing on your own. Then, instead of annotating, I think you should find already computed annotations of TEs.

The alternative to subscription-based repbase is dfam. And Just looking at their webpage they seem to have an annotation of all individual TEs in hg38. I bet my shoes there will be plenty of TE annotations in human in genome databases too.

---
Tags: repeat-elements, repeatmasker, transposable-elements
---


# Question
Title: Differentiating molecules based on peptide sequence? How to annotate?
I want to differentiate between classical class I and non classical class I MHC molecules in a model organism using well conserved structural features within classical MHC I molecules (eg intradomain disulfide bridges within class I a1 domain). To do this I need to align the peptide sequences of transcripts of that model organism with transcripts of classical MHC I human genes and see which genes from those transcripts have those structural features. So I am looking for a tool to align and annotate those structural features in a replicable manner

Additionally if you have any other bioinformatic ideas /computational analyses on how to find classical MHC Class I molecules on an organism which is poorly annotated that would highly appreciated!

# Answer
> 1 votes
Knowing almost nothing about the biology, I think that the most straightforward way to do this is:

1. train profile HMMs on your different classes of protein sequences (in human maybe, but ideally including MHCs from a large number of organisms)
2. translate your transcripts into proteins and run them against your trained profile HMMs
3. confirm that the profiles adequately discriminate the two classes and adjust if needed.

For an example of how to do this (on MHC), you can see here.

These approaches are then easily extensible to further organisms, as you have the profile HMMs.

It's possible that a simpler approach would consist of just finding hits against human MHC and then just checking for the existence of sequence patterns corresponding to those structural motifs, but I doubt that would work well in any significantly divergent organism.

---
Tags: sequence-alignment, protein-structure, features
---


title: How to create this function for lapply?

text:
# Question
Title: How to create this function for lapply?
I have two lists of dataframes. List `x` and list `my_list`. Each having four dataframes and their names are same for the two lists, i.e., `bray-curtis, Chao1, FB`, and `Shannon`.

The dataframes within the two lists have same element names, but they contain different information. This is because list `x` has been generated after some operations on list `my_list`.

Now I want to write a function and run that on list `x` with \`lapply. the command I want to write as a function is

```
forest(x$Chao1, slab = paste(my_list$Chao1$StudyID))

```

Here, the first `Chao1` is coming from list `x` while the second `Chao1` is coming from the list `my_list`. It works for a single dataframe but I am unable to write the codes with `lapply` so that it runs across all the dataframes of the lists.

The desired outputs are four plots from four dataframes of list `x`. Can anyone please help me?

Update: @gringer's approach works, but it saves all the plots in a single list. When I tried to save the images with following command it just saved the last image:

```
tiff(filename = "~/Desktop/image.tiff", width = 30, height = 30, units = "cm", res = 1200)

```

Is there a way to save all them individually?

Thanks

# Answer
> 2 votes
Assuming the lists are named, and you want to fetch the element with the same name from each list, I'm guessing something like this would work (or at least get you closer to a solution):

```
v <- lapply(names(x), function(var){
   forest(x[[var]], slab=(my_list[[var]])$StudyID)
})

```

When creating images that don't have a multi-page format (i.e. anything except PDF), and not changing the image name with each save, you need to add in a format specifier to tell R what to do with new created images. For example, this will add a three digit zero-padded number for each new TIFF image:

```
tiff(filename = "~/Desktop/image_%03d.tiff", width = 30, height = 30, units = "cm", res = 1200)

```

---
Tags: r
---


# Question
Title: What is the current state-of-the-art in assembling hybrid transcriptomes?
We are considering attempting de novo assembly of a species transcriptomes (i.e. without a reference genome) using the combined NGS outputs of Iso-seq and Illumina.

One example I saw (Li et al 2017), used the standard PacBio tools to assemble a transcriptome, followed by correcting the assembled sequence using a tool called `proovread`, followed by `cd-hit-est` and then `cogent`.

Would this be considered a reasonable pipeline or is there an alterative recommendation?

# Answer
> 1 votes
My recommendation would be to run *Trinity* using the `--long_reads` option, which allows you to provide an error-corrected fastq file for anchoring reads to a transcript isoform. The reads are then clustered per-isoform, and assembled, in a similar fashion to what's done with genome-guided Trinity:

https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running-Trinity

```
Trinity --seqType fq --max_memory 50G --long_reads corrected_pacbio_reads.fa --left reads_1.fq  --right reads_2.fq --CPU 6

```

*Disclaimer: I have contributed code to Trinity at some time in the distant past*

---
Tags: assembly, transcriptome, pacbio
---


# Question
Title: Pathifier bioconductor package recommends to use technical noise for an argument. How to determine it?
I'm running the pathifier approach against C2 pathway curated database for a specific microarray dataset.

As I was reading the documentation of the pathifier, in order to configure it properly for my dataset, I saw that for the `min_std` argument they suggest to use the technical noise.

> min\_std: The minimal allowed standard deviation of each gene. Genes with lower standard deviation are divided by min\_std instead of their actual standard deviation. (Recommended: set min\_std to be the technical noise).

But how am I going to calculate it ? At first run I just used the value `0.2254005` from the example they provide with that package, but this might not be the right for my dataset.

Update:

So from here it says: `The minimal standard error (min_std) was set as the first quartile of the standard deviation in the data`. How am I supposed to calculate this first quartile of the std? Any idea? Should I calculate the sd() for each gene and then find the first quartile?

What do you think about this function?

# Answer
*Answer from @llrs, converted from comment:*

Probably is a question for the maintainer but I would *guess* that is a variance of the gene in the control dataset. So, if in normal samples the gene has a standard deviation of 0.5 this is the expected too in the tumoral cells. calculate the sd only on the control samples. So if the healthy/control samples have such variation the altered ones should have these or more according to the article (if I understand it correctly).

`std <- apply(data[, controls], 1, sd, na.rm = TRUE); quantile(std, 0.25)` should be much faster and efficient

But I would check it in support.bioconductor with the authors/maintainers.
> 1 votes

---
Tags: statistics, microarray, pathway
---


# Question
Title: how do I predict bacterial small non-coding RNA for a specific mRNA?
I am working on Vibrio parahaemolyticus. I have a gene that it is possible might be regulated by small non-coding RNAs. How can I predict possible sRNAs that target the transcript of this gene. Most software offers the opposite of what I want, i.e., they predict targets for small RNAs!

I'm trying to predict bacterial small non-coding RNA in general (Support Vector Machine) but I got more than 8 thousand sequences within intergenic regions alone. Is there a code for aligning these sequences to my transcript?

# Answer
*Answer from @llrs, converted from comment:*

You could try to download RNAseq data from this bacteria and see if some of the RNA align to your sequence. You could try to blast your gene and look for the sequences that have low identity similarity (or classify them by type and then look what is each hit).
> 2 votes

---
Tags: bacteria
---


# Question
Title: How can I calculate loss of heterozygosity (LOH) in NGS sequencing data?
I'm analyzing a tumor sample and a healthy sample from the same patient. I want to use sciClone to look at tumor clonality. One input is the genomic regions that need to be excluded due to LOH. I have both WES and WGS.

Is there currently a "gold standard" for which tool will return the most robust regions?

The sciClone paper mentions VarScan2, APOLLOH, HMMScan? Are these still state of the art?

Is it feasible to calculate LOH programmatically without the need of external software?

# Answer
*Answer from @fabio-marroni, converted from comment:*

WARNING! Self-promotion! We published a freely-available tool for detection of LOH (and also redeuction fo homozygosity). Here you find the link to the paper and therein will find the link to the software. We compared with other tools that you might also try to see which one is the best performing. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5466433/
> 1 votes

---
Tags: cancer
---


# Question
Title: Getting exact mismatch data from multiple sequence alignment
TL;DR Can anyone point me in the direction of a tool that can be used to input different viral substrain reference sequences (fasta) as well as sequence data (fastq/bam/vcf) and find the most likely substrain of the sequenced data?

I am doing multiple sequence alignment of different viral subtype references and would like to get the exact positions and base each time there is a mismatch. I have tried using different tools. For example Clustal Omega can almost give me what I seek, but I have trouble getting the exact position and nucleotide for the mismatches (the full thing is ~8000bp long):

```
ACTACAATAATCCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAATTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG

*********** *********** ****** *****************************

```

E.g I would like something like:

```
Pos 12:
C
T
T
T
T
T
T
T
T
T
T

```

For each position with a mismatch.

A problem is that the sequences are slightly different lengths so something like a variant call format on the "main" strain might work, but I'm not sure how I can do this on multiple files of this length.

The end goal is actually to be able to call which substrain is present from NGS data. I imagine using variant calling on the "Main" strain, but getting from there to calling which substrain is in the data is a bit of a long way for me. If someone have any suggestions on tools for this, please enlighten me :)

Thank you very much!

# Answer
> 0 votes
I found a solution. For variant calling on multple sequence alignment I can recommend snp-sites. For anyone facing a similar problem of wanting to find substrains, here is an overview of useful tools: "Computational Methods for Strain-Level Microbial Detection in Colony and Metagenome Sequencing Data" (https://www.frontiersin.org/articles/10.3389/fmicb.2020.01925/full)

# Answer
> 0 votes
I have partly answered this for know! For anyone interested i use kalign to align all 10 strains to a fasta format, then I use snp-sites to output a vcf file that can show each mismatch in accordance to the first strain in the fasta. New problem is that the alignment "stretches" the length of this first strain, because some of the other strains are a bit longer. This means that the position of variants are with reference to the longest strain of the 10 and not the strain I want. It can work for now, but any solution to this, would be helpful.

---
Tags: multiple-sequence-alignment
---


title: Get Salmon mapping/alignment summary

text:
# Question
Title: Get Salmon mapping/alignment summary
With HISAT2, after the alignment of fastq files you get an alignment summary like this:

```
20000 reads; of these:
  20000 (100.00%) were unpaired; of these:
    1247 (6.24%) aligned 0 times
    18739 (93.69%) aligned exactly 1 time
    14 (0.07%) aligned >1 times
93.77% overall alignment rate

```

Is it possible to get this information with Salmon?

# Answer
> 6 votes
The `salmon` output log contains a line with `[info] Mapping rate` which allows to calculate how many reads went mapped or unmapped. Is that sufficient? It should be in the `logs` folder of the directory `salmon` creates per sample.

---
Tags: rna-seq, sequence-alignment, salmon
---


title: How do I install hammock?

text:
# Question
Title: How do I install hammock?
I'm working with phage peptides libraries. Currently I'm analysing a massive data set from phage selection, NGS readout, preprocessed already (fasta file). I want to analyse the sequences to cluster peptides according to their consensus motifs. I found that Hammock may be helpful to do this.

Here is the GitHub repo for Hammock: https://github.com/krejciadam/hammock. I downloaded the master, unpacked it, but I don't know how to use it.

I guess I need to compile, but I have no idea what files need to be compiled, how to do it, what software I need use to compile etc.

I work with Mac/terminal.

# Answer
> 2 votes
1. Download a release (Hammock\_v\_1.2.0.7z) from https://github.com/krejciadam/hammock/releases
2. Unzip
3. To run the first example:

```
cd /your_path/Hammock_v_1.2.0/examples/MUSI/
java -jar ../../dist/Hammock.jar full -i musi.fa

```

As mentioned in the comments, it's in the manual and any further issues need more information (errors etc).

# Answer
> 0 votes
*Answer from @ram-rs, converted from comment:*

There's a Hammock.jar file - there' no need to compile Hammock, although they recommend compiling dependencies. Have you tried running java -jar dist/Hammock.jar --help? Follow their steps from Example 1 - it should work.

---
Tags: clustering, software-installation
---


title: fastq file format unknown

text:
# Question
Title: fastq file format unknown
I have paired-end fastq files some of which seem to be in a weird format (from a collaborator, not a public database). When I `cat` the file I get what seem to me to be binary syntax but not sure?

However when I `cat` into the other pair-end the output seems fine

Any help/suggestion/solution would be greatly appreciated.

It does seem to be compressed with a missing suffix:

```
file pair1.fastq
   pair1.fastq: data.

```

# Answer
> 3 votes
This indeed looks like some sort of binary file. I have seen bunch of those when people renamed `.fastq.gz` to `.fastq` without actually uncompressing them. However, as you say, this is not your case, as `file pair1.fastq` would identify the file as `gzip compressed data`.

Looking at possible return values of `file` comments here, `data` basically mean the command could not figure out anything about it beside being some sort of binary data file, but I suspect it probably somehow broken file. Honestly, at this point I think you should check the source of the data. If it's from a public database, re-download the data. If it's your own sequencing, contact the sequencing provider and check the data with them. I don't see any way to rescue this file (which does not necerasrily mean there is nothing you can do).

---
Tags: sequence-alignment, ngs, trimming, bowtie2
---


# Question
Title: Calculating most abundant transcript from RNA-Seq data
`vcf2maf` uses `VEP` to annotate variants, and I believe selects the default Ensembl transcript to use for annotation. Sometimes the transcript that `VEP` selects is not the transcript I'm interested in, usually because the selected transcript is not the most highly expressed transcript in my tissue of interest (skin). `vcf2maf` allows you to provide a transcript override list so that VEP annotates the variant using the specified transcripts instead.

I have several skin samples sequenced with bulk RNA-Seq. I want to estimate the average abundance for each transcript across all samples and then use these abundances to rank transcripts from most to least abundant. Then I will use the most abundant transcript as the default VEP transcript. I plan to use `salmon` or `kallisto` to quantify transcript abundance. Should I use TPM or normalized counts to calculate average expression?

My initial thought is to use normalized counts (generated by DESeq2 from raw counts). Are there any problems with this approach? GTEx displays transcript abundance with average TPM, but I thought TPM was inappropriate to use across samples because it doesn't account for between sample differences.

Update: I forgot to mention I also tried using TPM ranks like @ATpoint describes. I haven't fully compared how this compares to transcripts identified by normalized counts, but the initial genes I checked showed good concordance between methods

# Answer
Transcript abundance quantification is a tricky topic since a read often could belong to several transcripts, so any "count" is a best guess as to which transcript it actually originates from. That being said, there are tools that can help you here:

1. `salmon` (as you mentioned) to quanitfy. Run it with `--numGibbsSamples 50` (or higher if your computer has the hardware. Also run with `-d` for the next tool
2. `terminus` is a tool that also originates from the COMBINE-lab, It's relatively experimental, but it applies well to your use case. It will take the salmon output and assign counts to transcript groups. If a set of transcripts can be clearly quantified above a certain confidence, that transcript will stay unique. If it is confounded with other transcripts, those will be placed in a transcript group and their counts will be summarised. More info here

You are right in that TPM, FPKM etc. are not appropriate measures to compare across samples. TMM, or DESeq2's median ratio method (what you get with `counts(dds, normalized=TRUE)` are a step in the right direction, but are still inappropriate across samples. I'd suggest this often posted blog post in the topic.

There are no perfect solutions here since RNA-Seq is inherently a relative quanitifcation method so this will always be an issue. IMO, use `DESeq2::varianceStabilizingTransformation()`, `DESeq2::rlog()` or `limma::voom()` after terminus and rank your transcripts that way.

Just keep the inherent transcript incertainties in mind.
> 5 votes

# Answer
If I read you correctly, you want to find the gene with the highest expression?

With RNA-seq data from multiple samples, I'd find this by ranking each gene by TPM/FPKM and taking the mean rank across samples. This method will tend to mitigate any large differences in a gene's expression between samples.

As you are comparing (RNA-seq) expression across genes you must take into account transcript length.
> 0 votes

---
Tags: rna-seq, deseq2, transcriptome, normalization, salmon
---


# Question
Title: Global alignment between two sequence X and Y with maximum number of identical matches
If 2 protein sequence X of length m and Y of length n and if there are several highest scoring global alignments, then I want to get one alignment that has largest number of columns in which a letter from X is aligned to an identical letter from Y. I need an algorithm to find alignment with maximum number of identical matches. Can anyone help me?

# Answer
There are a lot of possible ways to count the number of perfect matches in an alignment, and the details will depend on which alignment tool you are using.

But, if you use the built-in pairwise aligner in BioPython (`Bio.Align.PairwiseAligner()`, https://biopython.org/docs/1.75/api/Bio.Align.html), then the easiest way is to use the formatted string representation of the alignment, which uses the pipe/vertical-bar character ('|') to indicate perfect matches.

Example code:

```
from Bio import Align

aligner = Align.PairwiseAligner()
alignments = aligner.align("GAACT", "GAT")
max_matches = 0
for i, a in enumerate(alignments):
    align_str = f'{a}'
    align_str = align_str.split('\n')
    match_str = align_str[1]
    match_count = match_str.count('|')
    mismatch_count = match_str.count('.')
    print(f'Number of matches: {match_count}')
    print(f'Number of mismatches: {mismatch_count}')
    if match_count>max_matches:
        max_matches = match_count
        max_ind = i

best_alignment = alignments[max_ind]

```

In this example, both alignments have the same number of perfect patches (3), but that may not always be the case, depending on the input sequences and the scoring parameters used. The example just uses the default scoring parameters, but in general, you probably want to choose parameters appropriate for your type of alignment.
> 1 votes

---
Tags: sequence-alignment
---


# Question
Title: Glocal\semi-local\Hybride Globale-Local alignment with Python
I was looking for a simple way to do a glocal alignment. The case I have is I have a small sequence which should be find in a bigger one, thus typically a glocal alignment. Also I can not Install anything on the target machine (Python 3.6 with Biopython is already installed), all I can do is pip install.

My question is : are there any glocal alignment implementation out there ?

I was unable to found one and I don't feel proficient enough with Python yet to implement it myself.

# Answer
> 1 votes
As suggested in the comments, the built-in BioPython pairwise aligner should work for this. You can find a short sequence match in a longer sequence using a global alignment with the right score parameters (`query_end_... = 0`). The other scoring parameters are values that I've found to be good for aligning similar sequences.

Example code:

```
from Bio import Align

aligner = Align.PairwiseAligner()
aligner = Align.PairwiseAligner()
aligner.match_score = 5
aligner.mismatch_score = -9
aligner.mode = 'global'
aligner.target_internal_open_gap_score = -12
aligner.target_internal_extend_gap_score = -3
aligner.target_end_open_gap_score = -12
aligner.target_end_extend_gap_score = -3
aligner.query_internal_open_gap_score = -12
aligner.query_internal_extend_gap_score = -3

# Set the query_end penalties to zero to effectively implement a 'glocal' alignment as described in the question
# These are the default values, so you could skip the next two lines,
#    but this is relevant to the question, so I'm including them in the example.
aligner.query_end_open_gap_score = 0
aligner.query_end_extend_gap_score = 0

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATCCG")

for a in alignments:
    print(a)
print()
print()

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATACG")

for a in alignments:
    print(a)
print()
print()

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATACCG")

for a in alignments:
    print(a)

```

---
Tags: sequence-alignment, python
---


# Question
Title: How to calculate module-trait relationship when trait data is in binary format?
I have a dataset of 50 breast cancer samples. These samples are classified into four subtypes Lum A, Lum B, Her2 and Basal. I have been working with lncRNAs and protein-coding genes. To identify the functions of lncRNAs, I have used WGCNA through which I get the modules with protein-coding genes and lncRNAs. The expression data is stored in a matrix `datExpr`.

I have the trait information in `datTraits` like below:

I totally got 20 modules named with different colors stored in `dynamicColors`. Like below I first calculated eigenvalues and then used it and trait information to check the association between Modules and trait data.

```
MEs0 = moduleEigengenes(datExpr, dynamicColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

```

What is the optimal method to associate modules and trait using binary data?

1. Is Pearson correlation a right way to check the association?
2. Is there any other way to check the association between modules and trait (which represented like above in binary)? Details of methology would be welcome.

Patients are not genetically related and I have no information on environmental risk factors. I was thinking about logistic regression, but before doing that, I would like to take suggestions.

# Answer
> 1 votes
*Answer from @m, converted from comment:*

Personally I'd build a tree. Whats your question patients or the relationship between RNAseq 'signatures'? Its not the definitive analysis, but its good mining.

Pearson is doable but how many Pearson's? You'll end up with a 2x2 matrix. How will you resolve the matrix clustering? A maximum likelihood phylogeny solution would work here making use of the binary format. There are other approaches.

Spearman correlation would be better than Pearson's you end up with a 2x2 matrix.

*Answer from @peter-langfelder, converted from comment:*

You can use Pearson correlation, it is equivalent to carrying out an equal-variance t-test for the module eigengene between the two groups specified by the binary variable.

---
Tags: r, rna-seq, networks, wgcna, correlation
---