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

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.


Questions regarding Nanopore sequencing analysis I am new to Nanopore sequencing analysis. I have a couple of questions regarding it which are as follows:

How do I know if my fast5 file is multiread or single read file?
Is there a way to combine all the fast5 files into a single fast5 file?
How can I extract the flowcell type and sequencing kit information from the fast5 file? Is there a command that I can use to extract it? I know I can use the following command to retrieve it:

h5dump -a /UniqueGlobalKey/context_tags/sequencing_kit -a /UniqueGlobalKey/context_tags/flowcell_type <fast5 file>

However, it works well for single read file and not for multiread files. Please correct me if I am wrong. Thanks!

1. The quickest way to know if it's a single or multi fast5 is the file name. If it mentions channels and read numbers (e.g. ..._read_9_ch_471_..., then it's a single read file. Alternatively, the multi-read files have only read IDs at the top level (i.e. read_<32-character-string-with-additional-dashes>). Example:

$ h5ls <read_file_name>.fast5 | head
read_0115573a-042a-42a3-b4e2-f0338b4a6c66 Group
read_015b870d-8a03-4f22-a73e-fb9f833b94c7 Group
read_0167393a-dede-4a27-91af-494622b63261 Group
read_01e658de-853b-4f58-93ad-87fbd735254c Group
read_0202acc2-33f4-4d29-9dac-c7f80c758c86 Group
read_022d68f0-6f8e-4124-bfb0-11262efa4a67 Group
read_0258a3c9-4ffe-42ab-a42d-89466fcefbf5 Group
read_032da6db-e84a-43ec-8c25-bcaead992066 Group
read_03383754-be43-420b-87fc-60c2eb833963 Group
read_036b4824-bb0a-42dd-ade0-29d89112ff07 Group

2. ONT provide a single_to_multi program for converting to the multi-read format. Details can be found on this github page. Here's a quick example of how it could be run:

$ single_to_multi_fast5 --input_path fast5 --save_path fast5_multi \
    --filename_base multi_file_myRunID --batch_size 4000 --recursive

3a. For multi-fast5 files, the kit/cell attributes are in the context_tags group of the file (because a multi-fast5 file could contain fast5 elements from multiple different sequencing runs with different run parameters):

$ h5dump -g /read_<readID>/context_tags multi_file_myRunID_0.fast5 | \
  grep -A 11 -e 'sequencing_kit' -e 'flowcell_type'
   ATTRIBUTE "flowcell_type" {
      DATATYPE  H5T_STRING {
         STRSIZE 11;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "flo-min106"
      }
   }
--
   ATTRIBUTE "sequencing_kit" {
      DATATYPE  H5T_STRING {
         STRSIZE 11;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "sqk-pbk004"
      }
   }

3b. For single fast5 files (at least, in the most recent version of them), they are found in the UniqueGlobalKey group:

$ h5dump -g /UniqueGlobalKey single_file_read_X_ch_Y_strand.fast5 | \
  grep -A 11 -e 'sequencing_kit' -e 'flowcell_type'
      ATTRIBUTE "flowcell_type" {
         DATATYPE  H5T_STRING {
            STRSIZE 11;
            STRPAD H5T_STR_NULLTERM;
            CSET H5T_CSET_ASCII;
            CTYPE H5T_C_S1;
         }
         DATASPACE  SCALAR
         DATA {
         (0): "flo-min106"
         }
      }
--
      ATTRIBUTE "sequencing_kit" {
         DATATYPE  H5T_STRING {
            STRSIZE 11;
            STRPAD H5T_STR_NULLTERM;
            CSET H5T_CSET_ASCII;
            CTYPE H5T_C_S1;
         }
         DATASPACE  SCALAR
         DATA {
         (0): "sqk-pcs108"
         }
      }

[note that this file structure may change in the future]


Block wise protein imputation I am currently working on a dataset that contains 50 samples (10 samples * 5 blocks). The features of the date set are:

The data is perfectly balanced between blocks, with equal treatment representation in each block. Each block contains 2 control samples (CTL) that are pools of all other treatment samples. These are technical replicates that are identical across multiplexes.

There is a lot of missing data (~40%),

most of the missing data is missing by block, meaning that the protein is present in some blocks, but not others.

When looking at the expression of these missing proteins in the other blocks, they tend to be around the median value for the experiment, which suggests to me that these values are missing completely at random.

I could proceed using complete cases, but then I'm missing out on a large portion of data. I am considering using an imputation method, but I'm concerned that (due to the nature of the missing data) that I will just be introducing a block effect into the data.
Any input on potential methods for imputation would be greatly appreciated.
Experimental design:
SampleNumber    Multiplex   Treatment
1   A   CTL
2   A   CTL
3   A   TMT1
4   A   TMT1
5   A   TMT2
6   A   TMT2
7   A   TMT3
8   A   TMT3
9   A   TMT4
10  A   TMT4
11  B   CTL
12  B   CTL
13  B   TMT1
14  B   TMT1
15  B   TMT2
16  B   TMT2
17  B   TMT3
18  B   TMT3
19  B   TMT4
20  B   TMT4
21  C   CTL
22  C   CTL
23  C   TMT1
24  C   TMT1
25  C   TMT2
26  C   TMT2
27  C   TMT3
28  C   TMT3
29  C   TMT4
30  C   TMT4
31  D   CTL
32  D   CTL
33  D   TMT1
34  D   TMT1
35  D   TMT2
36  D   TMT2
37  D   TMT3
38  D   TMT3
39  D   TMT4
40  D   TMT4
41  E   CTL
42  E   CTL
43  E   TMT1
44  E   TMT1
45  E   TMT2
46  E   TMT2
47  E   TMT3
48  E   TMT3
49  E   TMT4
50  E   TMT4

The data is in a matrix format, columns = samples and rows = individual proteins, similar to gene expression.
The controls are technical replicates, but they are pools of all other 40 samples. This is for TMT proteomics, where the pools are used by the MS to identify peptides across the multiplexes, not necessarily for multiplex correction. I've included multiplex (as well as treatment) in the linear model to identify changes associated with the treatment, while removing the effect of the multiplex. I have confirmed (with PCA and PVCA) that the multiplex effect is negligible after including it in the lm. My problem at the moment is the missing data.
I calculated the missing percentage from the data matrix by taking the number of missing values divided by the total number of potential values. he data is in a matrix format, columns = samples and rows = individual proteins, similar to gene expression. To calculate the missing values, I divided the total number of NA's by the total number of values in the matrix (number of samples * number of proteins). Missing data is common in proteomics, but normally this can be attributed to something, most often being close to the limit of detection. But that is not the case here, if a value is missing, it is missing in the whole multiplex, including the pool, where we know it should be found.
The missing values in this case are random, and this is definitely due to the MS. It will select the first ~65k peptides that it identifies to trigger an additional MS event, so if it doesn't see that first peptide, it will be missing from the entire multiplex. There are methods for imputation of data that is missing at random, but I'm concerned that it will just be imputing by multiplex.

40% missing data is huge. Missing data analysis is complicated on the underlying distribution. If the data set is periodic then missing data periodicity is needed. Non-periodic data can be solved using regression, I think SciPy has an automated method for this, but the volume of missing data makes this approach complex. If you can bring the missing data down to 20% by compartmenatlising the data set that would be the starting point.


How should the SAM MD tag match the CIGAR string? I am trying to understand how the MD:Z tag is used. The following is from the SAM Optional Fields Specification, which gives an example but is not thorough.

  The MD field aims to achieve SNP/indel calling without looking at the
  reference. For example, a string ‘10A5^AC6’ means from the leftmost
  reference base in the alignment, there are 10 matches followed by an A
  on the reference which is different from the aligned read base; the
  next 5 reference bases are matches followed by a 2bp deletion from the
  reference; the deleted sequence is AC; the last 6 bases are matches.
  The MD field ought to match the CIGAR string.

Suppose I have a read that I want to soft clip at both ends. If the read starts with CIGAR 100M and MD 50G49, and I want to change the CIGAR to 7S86M7S, what should the MD field become?

Is there a more complete explanation of what appears in the MD tag field?

The MD string doesn't apply to soft or hard clipped regions,. so your example read becomes 43G42. Since variant calling using something simplistic like this is only ever going to use the aligned portion, I guess there was never much reason to bother with adding in soft-clipping information. This also means that you need to parse the CIGAR string if you're using the MD string.


What should be the numbers in the fastSTRUCTURE file? The structure file for fastSTRUCTURE seem to be a table with rows and columns but the cells are all numbers, while my genotypes are simply bases. How could I know what I should insert there?

Numbers, usually in the form of 0/1/2, is a common way many formats represent bases.
For example, if a SNPs has bases which are A/T, then these are encoded as 0/1, with the reference allele/base taking the 0 and the alternate allele/base taking 1.
However, if you read the documentation, it says:

Specifically, rows in the data file correspond to samples, with two
rows per sample (note that only diploids are handled by this
software), and columns correspond to SNPs. The first 6 columns of the
file will be ignored; these typically would include IDs, metadata,
etc. This software only handles bi-allelic loci. The two alleles at
each locus can be encoded as desired; however, missing data should be
encoded as `-9'.

So you should be able to use your base letters here and it will work OK as long as the columns and rows are formatted correctly.


Tools for comparing/visualizing FASTAs? I have two FASTAs/assemblies of the same species, but the bases are somewhat different. I would like to explore this.

What tools/methods exist to compare two FASTAs and see the difference in structure between them? Are there visualization tools available, something like IGV?

I think users would need to do a multiple-sequence alignment between these FASTAs, and then check for the differences between them. However, I don't know if a "standard approach" exists to do this.

You can try Mauve multiple genome aligner.


How to find all variable-length seqs with an exact 5' and 3' match in a FASTA file Context

I am interested in finding all of the promotors specific to a particular sigma factor. I have identified the -35 and -10 sites from the literature, bold denotes -10, -35, binding sites:

  -35 site               -10 site
5' TTTACAtatttatttcagacaacGTCTTT 3'
   ^^^^^^                 ^^^^^^

The CAPITALIZED and ^emphasized^ 5' and 3' sequence components are the conserved sequence that must have an exact match. The lower-case bases are the variable 10-25bp region for which the sequence does not matter.

The sigma factor is from C. difficile and I am particularly interested in searching C. acetobutylicum, however, it would also be interesting to see where else this sequence crops up. Obviously such a short sequence will arise often by random chance, so it is important to define the spacer region, which might be anywhere between 15-25 nucleotides and may be any sequence.

In particular, I would like to find these binding regions from this sequence and this one.

Desired output

A text file that contains the locations and sequences of all of the variable-width sequences as defined above. Something like:

seq   start  stop  gap  dir  seq
seq1  1      22    10   F    TTTACAagtcagatctGTCTTT
seq2  2      23    11   F    TTTACAagtcagactctGTCTTT
seq2  122    101   10   R    AAAGACatgctagccaTGTAAA

Thank you.

python 3 solution

This python program will get you what I think you want. It is not elegant given that it (A) stores the whole sequence in memory and (B) spends most of its time slicing strings rather than working with a buffer. This program is not recommended for FASTA files for which the largest sequence will not fit well in RAM.

All that you need to do to modify it for other genomes or other binding sites is to change the three lines that begin with forseq, revseq, and filename.

from Bio import SeqIO
from Bio.Seq import Seq

forseq = Seq("tttaca".upper())
revseq = Seq("gtcttt".upper())
filename = "NC_003030.fasta"

forcom = forseq.reverse_complement()
revcom = revseq.reverse_complement()
print("seq\tstart\tstop\tgap\tdir\tseq")
for record in SeqIO.parse(filename, "fasta"):
    for i in range(len(record.seq)):
        for v1, v2, direc in [ (forseq, revseq, "F"), (revcom, forcom, "R") ]:
            if record.seq[i:i+len(v1)] == v1:
                for j in range(10,25):
                    if record.seq[i+len(v1)+j:i+len(v2)+j+len(v2)] == v2:
                        print("{}\t{}\t{}\t{}\t{}\t{}".format(
                            record.id, i,
                            i + len(v1)+j+len(v2)-1,
                            j, direc,
                            record.seq[i: i+len(v1) + j + len(v2)] ))

Here are the first few lines of example output from the first fasta file that you linked.

seq          start    stop     gap  dir  seq
NC_003030.1  97512    97543    20   F    TTTACAAAAGAAAAGATATGCATAATGTCTTT
NC_003030.1  283173   283205   21   F    TTTACACTTATGTTTATATAGTTCAACGTCTTT
NC_003030.1  303104   303125   10   R    AAAGACGCAAGCATAATGTAAA
NC_003030.1  1127691  1127725  23   F    TTTACAAGGTGTTACTGATAGAACAAAGGGTCTTT
NC_003030.1  1301958  1301988  19   F    TTTACATAAAAAGATGATAGTTTTTGTCTTT

python 3 with redundancy

If you need sequence redundancy, this version works. I used the redundancy codes available from www.reverse-complement.com You may use the redundant characters in [ R, Y, S, W, K, M, B, V, D, H, N ] or lowercase variants. I'm not sure how to handle the '-' character at the moment but am open to suggestions. Code and results below.

Again, the algorithm isn't great given all the for loops so other answers are appreciated!

from Bio.Seq import Seq

forseq = "tttaca" # Just use python strings here. No Seq() as above
revseq = "ntnntn"
filename = "NC_003030.fasta"

revcomp = {'A':'T', 'C':'G', 'G':'C', 'T':'A',
           'R':'Y', 'Y':'R', 'S':'S', 'W':'W',
           'K':'M', 'M':'K', 'B':'V', 'V':'B',
           'D':'H', 'H':'D', 'N':'N'}
redund = {'A':['A'], 'C':['C'], 'G':['G'], 'T':['T'],
          'R':['A','G'], 'Y':['C','T'],
          'S':['C','G'], 'W':['A','T'],
          'K':['G','T'], 'M':['A','C'],
          'B':['C','G','T'], 'V':['A','C','G'],
          'D':['A','G','T'], 'H':['A','C','T'],
          'N':['A','C','G','T']}

def revcom(seq):
    return "".join(revcomp.get(base, base) for base in reversed(seq))

def seqmatch(query, reference):
    """Checks if the query sequence (perhaps a promoter with redundant bases)
    matches the reference. Both query and reference must be the same length.

    Returns true if there is a match (redundancy allowed), returns false
    otherwise.
    """
    assert len(query) == len(reference)
    for i in range(len(query)):
        if reference[i] in redund[query[i]]:
            pass
        else:
            return False
    return True

forseq = forseq.upper()
revseq = revseq.upper()
forcom = revcom(forseq)
revcom = revcom(revseq)
print("seq\tstart\tstop\tgap\tdir\tseq")
for record in SeqIO.parse(filename, "fasta"):
    for i in range(len(record.seq)):
        for v1, v2, direc in [ (forseq, revseq, "F"), (revcom, forcom, "R") ]:
            if seqmatch( v1, record.seq[i:i+len(v1)] ):
                for j in range(10,25):
                    if seqmatch( v2, record.seq[i+len(v1)+j:i+len(v2)+j+len(v2)] ):
                        print("{}\t{}\t{}\t{}\t{}\t{}".format(
                            record.id, i,
                            i + len(v1)+j+len(v2)-1,
                            j, direc,
                            record.seq[i: i+len(v1) + j + len(v2)] ))

Here is some sample output:

seq          start  stop  gap  dir  seq
NC_003030.1  733    767   23   R    GATGAAGATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  736    767   20   R    GAAGATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  739    767   17   R    GATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  742    767   14   R    CAAGAAACCGATACAAACAATGTAAA
NC_003030.1  743    767   13   R    AAGAAACCGATACAAACAATGTAAA
NC_003030.1  2665   2697  21   R    TATTAAAATTTTAGAGGATGATGGTGATGTAAA
NC_003030.1  2668   2697  18   R    TAAAATTTTAGAGGATGATGGTGATGTAAA
NC_003030.1  3798   3822  13   R    CACCAGAAGACTTAAAAATTGTAAA
NC_003030.1  3801   3822  10   R    CAGAAGACTTAAAAATTGTAAA


Mixcr batch processing I have a folder with 150 pairs of fastq files from an illumina sequence run. I need to precess them with mixcr, how can I do this in a batch with a single command?

Here's how I'd approach this:

for READ1 in /Documents/Data/Input/*R1*
do
    READ2=${READ1/R1/R2}
    mixcr analyze amplicon --s hs \
        --starting-material dna \
        --5-end v-primers \
        --3-end j-primers \
        --adapters no-adapters \
        --receptor-type trb \
        --only-productive "$READ1" "$READ2"
done

Basically, loop over all the R1 FASTQ files in a directory, and calculate the R2 names as needed.

This is using BASH string manipulation, and assumes that substituting R2 for R1 will correctly calculate the R2 filename from R1. It also assumes that all the FASTQ files are in a single directory.


Chunk alignment in a name sorted bam for parallel processing I have a bam file with 1 billion alignment reads of which there are 700 million unique reads. I want to split the alignments into chunks for parallel-processing. Multi-alignments of the same read should be in the same chunk. I want to use pysam to solve the problem and would be very grateful for suggestions.

python will be to slow for this job. Here's a awk solution. One need to sort by read name and take track over the number of reads per chunk. If the number is reached the next chunk can only created of the read name is different.

$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
  BEGIN { part=0; line=n }
  /^@/ {header = header$0"\n"; next;}
  { if( line>=n && $1!=last_read ) {part++; line=1; printf header > part".sam"; print $0 >> part".sam" }
    else { print $0 >> part".sam"; line++; }
    last_read = $1;
  }'

The output will be in sam format.

You can controll the number of reads per chunk by setting the n value for the awk command.

EDIT:

A second version that will output in bam:

$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
BEGIN { part=0; line=n }
/^@/ {header = header$0"\n"; next;}
{ if( line>=n && $1!=last_read ) {print part,line; part++; line=1;}
  print line==1 ? header""$0 : $0 | "samtools view -b -o "part".bam"
  last_read = $1;
  line++;
}


Why do molecular generation models maximize “penalized logP” as a measure of drug-likeliness? I found that Lipinski's rule of five states that Log P (octanol-water partition coefficient, lipophilicity measure) usually should not exceed 5.
Many papers about drug discovery machine learning models tell about maximization of "penalized logP", but the paper they refer in the end does not contain any information on it. Let me show that.
Let's check 1, here we read about penalized log P:

The second task is to produce novel molecules with desired properties. Following (Kusner et al., 2017), our target chemical
property y(·)is octanol-water partition coefficients (logP) penalized
by the synthetic accessibility (SA)score and number of long cycles
y(m) = logP(m)−SA(m)−cycle(m)

where cycle(m)counts the number of rings that have more than six atoms.

Then we go into Kusner et al., 2017:

For the second optimization problem, we follow (Gómez-Bombarelli et al., 2016b) and optimize the drug properties of molecules. Our goal
is to maximize the water-octanol partition coefficient (logP), an
important metric in drug design that characterizes the drug-likeness
of a molecule.As in Gómez-Bombarelli et al. (2016b) we consider a
penalized logP score that takes into account other molecular
properties such as ring size and synthetic accessibility (Ertl &
Schuffenhauer, 2009).

Here is Gómez-Bombarelli, 2016 and it seems to contain nothing on the topic (I can be mistaken, please point me then), except

The objective we chose to optimize was 5×QED−SAS, where QED is the Quantitative Estimation of Drug-likeness (QED),37 and SAS is the
Synthetic Accessibility score.36Thisobjective represents a rough
estimate of finding the most drug-like molecule that is also easy to
synthesize.

Other examples of the papers also leading to the same papers are 2, 3, etc.
What does penalized logP mean for drug-likeliness and why it should be maximized?

A positive logP is good (hydrophobic), but high results in ADMET problems (cytochromes break overly hydrophobic compounds) and aggregation (gunk, tar, crud are words happily thrown around). However, in most cases the logP is not greater than +3 because good binding is result of hydrogen bonds etc —some bad algorithms will suggest a polyaromatic hydrocarbon that makes lots of pi-pi interactions but these are few.
Confusingly, a high synthetic accessibility score is bad (hard). 3 is the average and 10 is basically something impossible. This score is how hard it is to make (i.e. pricier), which is one of the major drivers.
The number of cycles is number of rings. Even if you get pi-pi, cation-pi and sulfur-pi bonds you also get aggregation due to stacking (a crud). So is somewhat bad.
So basically that score is
high=good minus high=bad minus high=bad

If difference in relative Gibbs energy of binding is thrown it, negative ∆∆G is good, so that would be subtracted. As its atom count dependent ligand efficiency is better, but that ranges from 0 to -1.2 kcal/mol. So would need scaling.


Nextflow: Can I parameterize similar processes inside a for loop? I have set of ~5 processes which perform very similar tasks such that their differences can be parameterized. Rather than putting 5 items through some channel (which seems complicated), or defining 5 separate processes, it appears to be possible in Nextflow to declare the processes in a for loop, and output into an array of channels.
Below is a prototypical example and it seems to work as intended. (I understand this could be done differently with some channel cross product operator but I am just giving this as an example of a process in a for loop.)
Is there anything wrong with doing this? Is there a better way to do it?
Bonus question - the name of the process below is "foo". I would like to have the index reflected in the name i.e. "foo_$i" but this doesn't seem possible by default - any way to do this??
some_inputs = ['prot','dna', 'rna']

some_parameters = [5,3,2,7,1]

// Set up an array of channels to output to
receivers = new ArrayList<Channel>()
for (i = 0; i < 5; i++) {
    receivers.add(Channel.create())
}

// Iterate 5 times to create 5 processes each with different parameters
for (i = 0; i < 5; i++) {

    // WOW - process can be defined inside a for loop!!
    process "foo" {

      // Extract parameters here so correct $i value is used
      p = some_parameters[i]

      input:
      val x from Channel.from(some_inputs)

      output:
      val "$x $p" into receivers[i] // Note that correct index $i seems to be used here

      """
      """
    }
}

// Confirm that each of the processes outputted independently with the correct parameter
receivers[0].view { "Received[0]: $it" }
receivers[1].view { "Received[1]: $it" }
receivers[2].view { "Received[2]: $it" }
receivers[3].view { "Received[3]: $it" }
receivers[4].view { "Received[4]: $it" }

Output - note that the parameter is consistent within each input channel i.e. working correctly:
N E X T F L O W  ~  version 20.10.0
Launching `play.nf` [fabulous_keller] - revision: c89e2b2281
executor >  local (15)
[10/b477ec] process > foo (3) [100%] 3 of 3 ✔
[16/8282ba] process > foo (2) [100%] 3 of 3 ✔
[db/f7327e] process > foo (1) [100%] 3 of 3 ✔
[fb/886f6f] process > foo (1) [100%] 3 of 3 ✔
[46/c34ed3] process > foo (1) [100%] 3 of 3 ✔
Received[1]: rna 3
Received[1]: prot 3
Received[4]: rna 1
Received[2]: rna 2
Received[0]: dna 5
Received[2]: dna 2
Received[2]: prot 2
Received[0]: prot 5
Received[3]: rna 7
Received[0]: rna 5
Received[1]: dna 3
Received[3]: dna 7
Received[4]: dna 1
Received[3]: prot 7
Received[4]: prot 1

Use channels, this is what nextflow's entire design is based on. I'm honestly surprised that spawning 5 processes with the same name in a loop even works.
Turn your parameters into regular channels, use the combine operator (https://www.nextflow.io/docs/latest/operator.html#combine) to transform the 2 parameter channels into a channel that has each combination of params/inputs. Feed this channel into the  foo  process. Afterward use branch (https://www.nextflow.io/docs/latest/operator.html#branch) to split your foo output channel into 5 different channels.
typechan = Channel.from(5,3,2,7,1)
numchan = Channel.from('prot','dna', 'rna')

cartesian = numchan.combine(typechan)

process foo {
  input:
  val numtype from cartesian
  output:
  val numtype into bar
  script:
  """
  echo $numtype
  """
}

bar.branch {
  one: it[1] == 1
  two: it[1] == 2
  three: it[1] == 3
  five: it[1] == 5
  seven : it[1] == 7
}.set {result}
result.one.view()
```


Best distance parameter for estimating physical interaction between residues in a PDB file We can calculate the distance between residues in a PDB file regarding different parameters like closest atoms, alpha carbon, beta carbon, centroid and etc. Which one of these parameters are better to show physical interaction between residues in a PDB file?

A popular (and I would say, respected and trusted) website is PDBsum https://www.ebi.ac.uk/pdbsum (which also has a Wikipedia article about it: https://en.wikipedia.org/wiki/PDBsum)

They measure protein-protein (prot-prot) contacts as any Nitrogen, Carbon, or Oxygen element atoms (N, C, O) within 4.0 Angstroms (Å) when measuring the 3d euclidean distance.

If your atomic contacts are not close enough, the amino acids may not actually interact at all, even though they theoretically could.  Since in PPIs you generally want to use reliable ground-truth information, it is better to use cautious parameters.  Often, in PPIs, the answer to "does it interact?" is not "yes" or "no" but rather "yes" or "don't know / unknown".


In your ATAC data frame, you have chr1   16103  16348. But in your merged data frame, you have chr1   16104  16348. Is it possible that you a is 1-based whereas your ATAC data frame is 0-based? This can happen in bioinformatics and can be annoying. This might be why you have a lot of NA in your merged data frame. You should subtract 1 from the start position of your a data frame while keeping the end position unchanged. Then try to merge again.

You should always be aware of your data format to know whether it is 1-based or 0-based.


How can I compare two bed files? I am trying to translate (lift over) bed files describing genomic regions from hg37 to hg38. I have tried both UCSC's LiftOver tool and CrossMap but saw that they give me different results. I therefore need a way of assessing how correct the results are.

I decided to test by downloading a bed file of the RefSeq genes for hg37 and one for hg38 from UCSC. Then, I run liftOver to translate the hg37 file to hg39 coordinates and now I want to compare the results of the liftover to the bed file I downloaded for hg38.

How can I usefully compare these files? I am hoping for a tool that can report how similar the two are. Ideally, I would like to see three things in the output:

Matching regions: those that are shown with identical positions in the liftover results and the UCSC bed file.
Non-matching regions: those whose coordinates do not match between the two files.
(if possible) fuzzy matches: those that are not identical but pretty close (off by a few nucleotides).

I thought it would be easy to do with a little awk script but it gets quite complicated because certain RefSeq genes are mapped to multiple positions (e.g. MIR4454). And the 3rd requirement (fuzzy) is not trivial to do in awk. I can live without that though, if a tool can give me the other two.

From a set operations viewpoint, consider your hg37 regions as a "reference" set and hg38 as a "map" set.

If you have these sets in BED format, BEDOPS lets you do set operations on them with various overlap criteria.

Specifically, you can use BEDOPS bedmap to map hg38 regions to hg37 regions. The bedmap tool lets you assign overlap criteria from as little as one base to as much as full, mutual overlap.

Once you do mapping, you can do operations, like counting the number of regions that meet the overlap criteria between reference and map sets, and calculate percentages.

For example, use --exact to get exactly-matching regions between hg37 and hg38:

$ bedmap --count --exact hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > exact_counts.txt

You can do fractional ("fuzzy") overlaps, where (for example) at least 50% of an hg38 region has to overlap an hg37 region:

$ bedmap --count --fraction-map 0.5 hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > fuzzy_0p5_counts.txt

Getting a percentage out of a "fuzzy" result can be complicated by overlaps between hg38 regions that overlap a reference hg37 region. A more complicated procedure can be implemented with the --echo-map option and some command-line work to implement logic to deal with that case more stringently.

To get a straight-up count of no overlaps, build the inverse — or mask — of your hg38 regions, using Kent tools fetchChromSizes and BEDOPS bedops --difference:

$ fetchChromSizes hg38 | awk '{print $1"\t0\t"$2}' | sort-bed - > hg38.bounds.bed
$ bedops --difference hg38.bounds.bed hg38.bed > hg38.masked.bed

The regions in hg38.masked.bed make up all of the genomic space in hg38, except for your lift-over regions.

Apply --count --fraction-ref 1.0 to count the number of reference (hg37) elements which completely overlap these masked hg38 regions:

$ bedmap --count --fraction-ref 1.0 hg37.bed hg38.masked.bed | awk '{s+=$1}END{print s;}' > no_counts.txt

The documentation for bedmap describes these overlap criteria and operations in more detail.


Coronavirus RNA structures? Is there anything known about the RNA structures of coronaviruses? More specifically - do they have any interesting known structures in the translatable region, like RRE of HIV or the double loops in flaviviruses?

Update
Here is a recent development on the side of structure prediction.

To add a more complete answer: the current coronavirus is closely related to the SARS virus that caused the outbreak in 2004, and on which much research has been done.

Here is a general review of the coronavirus epidemiology, life cycle etc.

I haven't found yet any materials about the RNA structures in the translatable region, however the structures in the 3' and 5' terminal regions are well studied and are known to be crucial for initiating viral replication and for its packaging.

This recent article addresses in details the RNA structures in 5' terminal of human coronaviruses, whereas this older article describes the structures both in 3' and 5' terminals.

Finally, the interesting algorithmic developments for the link that I added to the question are described here.


Interpreting contig alignments to a reference genome I have applied two de novo genome assembly tools to data from the run SRR12707453, corresponding to a phage (I downloaded the data and have no relation with the authors of the study).
Using rnaSPAdes I get many long contigs (>50.000 bp). When analyzing the contigs using Quast against a reference genome (MN813694.1) most long contigs are discarded and the seventh-longest one is the only one selected. If I blast the other contigs I find correspondence with other phages, which makes me conclude that many phages were present in the same probe.
This seventh contig has 47686 bp length, while the reference has 47713 bp length. But they are not equivalently arranged. Apparently, the genome was circularized and then it was broken around the position 32400 in the reference. Also, the sequence does not match between positions 45000 and 45951 (in the reference), and a little misassembly seems to occur around position 46300.
See here some screen captures for the Icarus visualization obtained with Quast:

In contrast, using Quast with the contigs obtained from the other assembly tool, I get some differences. There is also a contig with almost the same size as the reference genome (47573 bp) that has exactly the same structure than the contig obtained with rnaSPAdes: circularized, broken around position 32400, does not match between position 45000 and 45951, and small misassemblie aronund position 46300. But there are also some contigs in which there is no break around position 32400.

May I interpret that these contigs indicate that the phage genome is sometimes circularized and sometimes broken at position ~32400 and that's why some contigs don't find a sequence break at that position?
Am I interpreting anything wrong? Any other conclusion can be obtained from the data?
You can (temporaly) download here the Quast results.

some additional information to your request.
https://ftp.ncbi.nlm.nih.gov/genomes/Viruses/FamilyPhylogeneticTree/Graph/SRR12707453.MN813694.1.gfa
it contain assembled by reference graph
using ncbi tool saute(skesa) which will be public later.
install bandage app and load this graph
https://rrwick.github.io/Bandage/
picture looks intriguing,I would say , run contain many errors located in one position.
Reason? I don’t really know.


Coordinates to Genome R Studio Script? I have a very, very long list of genome coordinates of areas. I need to gather a sequence in order to compare it to another database of sequences that I have. How do I produce the genome sequences from these coordinates in R?

I am very new to this and need a detailed answer.

This doesn't use R, but a simple and quite convenient method is to use bedtools getfasta with your fasta file. You'll need to make a BED file out of the RepeatMasker output:

awk 'BEGIN{OFS="\t"}{name=sprintf("%s_%s", $10, $11); strand="+"; if($9 == "C"){strand="-"}; print $5,$6-1,$7,name,".", strand}' repeatmasker.output > rmsk.bed

Then, given that and your genome.fa fasta:

bedtools getfasta -name -s -fi genome.fa -bed rmsk.bed -fo output.fa

This will produce a multiline fasta file names like AluSx_SINE/Alu (due to the -name option) where the sequences are reverse complemented if appropriate.

You can install bedtools either via conda (conda install -c bioconda bedtools) or github.


What is the difference between a transcriptome and a genome? I have a computer engineering background, not biology.

I started working on a bioinformatics project recently, which involves de-novo assembly. I came to know the terms Transcriptome and Genome, but I cannot identify the difference between these two.

I know a transcriptome is the set of all messenger RNA molecules in a cell, but am not sure how this is different from a genome.

In brief, the  “genome”  is the collection of all  DNA  present  in  the  nucleus  and  the  mitochondria of a  somatic  cell. The initial product of genome expression is the “transcriptome”, a collection of RNA molecules derived from those genes.


When running ChIPQC bioconductor package after peak calling, should I use sorted BAM files or unsorted BAM files? I have called peaks using MACS2 and now I'm trying to run ChIPQC bioconductor package. Can someone please tell me should I use sorted BAM files or unsorted BAM files? I have created sorted BAM files and from that indexes for BAM files. When I was trying run ChIPQC earlier, I got an error. Then I figured out that is because I don't have index files. Then I created indexes after sorting BAM files.

You need to provide sorted bam files, ideally with an index. If you run ChIPQCsample, it calls ChIPQC:::sampleQC which iterates through the chromosomes or contigs, using the index, and the absence of the index will crash it.
We can use an example, first in linux:
wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibTfbs/wgEncodeHaibTfbsA549RxlchV0416101Etoh02AlnRep1.bam -O Etoh02AlnRep1.bam

samtools index ./Etoh02AlnRep1.bam
# s
samtools view -uh Etoh02AlnRep1.bam chr1 > sorted.bam
samtools index sorted.bam
# quick way to unsort
samtools sort -n -o unsorted.bam sorted.bam
# it fails to index
samtools index unsorted.bam
[E::hts_idx_push] Unsorted positions on sequence #1: 246814554 followed by 42922635
samtools index: failed to create index for "unsorted.bam"

Then in R:
library(ChIPQC)
ChIPQCsample("unsorted.bam")
Bam file has 25 contigs
Error in value[[3L]](cond) :
  failed to open BamFile: failed to load BAM index
  file: unsorted.bam

ChIPQCsample("sorted.bam",annotation=NULL)
Bam file has 25 contigs
Calculating coverage histogram for chr1

Calculating SSD for chr1

Calculating unique positions per strand for chr1

Calculating shift for chr1

 300 / 300
Counting reads in features for chr1

....


Publicly available, free, complete database for antibiotics names and classes? This is a tough one I think: is there a publicly available, up-to-date, free, complete database for antibiotics names and classes?

I am specifically looking for information like, e.g., cefpirome (is a) cephalosporin.

I've looked in a couple of different places, but nothing really comes close. E.g.:

CARD comes closest, but is not free (companies need to pay a license) and not complete/up-to-date (e.g., cefpirome is missing)
@Karel mentioned KEGG in his answer, but this also needs a license (companies for sure, academics if they download data via FTP).
Novel Antiobiotics Database has no class ... and stops in 2003
The Antimicrobial Index is web only, and terms of use apparently prohibit any kind of web-scraping.
A suggestion by @Pierre on Twitter was to try Antimicrobial stubs on Wikipedia. Good idea, but the class is not stored in the ChemBox and must be free-text parsed. Or writing a script which evaluates the Wikipedia category tags.
I drew blank at different government agencies (EMA, FDA, EFSA, NIH, etc.). The closest I found was from the WHO a list of critically important antibiotics. But not complete, and in a PDF.

Any pointer appreciated.

ChEBI

It should be possible to get this information from the ChEBI database, see the exported tables. You could download the ontologies (in OWL / OBO), parse them using some ad-hoc parser or using a dedicated library (e.g., Pronto), build a directed acyclic graph based on the is_a relationship, and extract the antibiotic subtree.

KEGG

You could also use KEGG: Kyoto Encyclopedia of Genes and Genomes, but a subscription is probably required. The exported drug table should be easy to parse (e.g., in Python).

For a quick verification of the completeness, you can download the table, filter out everything except antibacterial agents and check the number of records:

curl "ftp://ftp.genome.jp/pub/kegg/medicus/drug/drug" > keg_all.txt

cat keg_all.txt \
   | perl -pe 's/\n/BREAK/g' \
   | perl -pe 's@///BREAK@///\n@g' \
   | grep 'CLASS       Antibacterial' \
   | perl -pe 's/BREAK/\n/g' \
   > keg_antibacterial.txt

cat keg_antibacterial.txt | grep "///" | wc -l

There are 646 records. However, some of them are almost identical (e.g., only 1 atom difference).

If you add | grep -i "antibiotic" \ to the pipeline, then you get only 445 records.


Low Fraction of usable antibody reads in CiteSeq we performed a combined gene expression and CiteSeq experiment with the 10x VDJ kit and 20 conjugated antibodies and sequenced on hiseq. I used cellranger to process the sequencing output. The cellranger summary shows overall good values except for a low percentage of "Antibody Reads Usable" and a low fraction of "Antibody Reads in Cells". How is it possible that so few antibody reads are usable? Does this mean that most antibodies were not bound to the cell surface?

We have also seen a very similar issue and are also investigating why. Possibly of note, I tried forcing the cellranger program to assign a very large number of cells (10x the actual number of cells we had) using the --force cells option. This appeared to increase the % of usable Ab reads and Ab reads per cell metrics, suggesting that these reads may be associated with barcodes that have relatively few UMIs?


How are snakemake's --cluster and --drmaa options implemented? I'm fairly new to snakemake and I'm trying to understand the difference between the --cluster and --drmaa flags, both of which are used to submit jobs to compute clusters/nodes.

The docs give a few hints about the advantages of using --drmaa here:

  If your cluster system supports DRMAA, Snakemake can make use of that
  to increase the control over jobs. E.g. jobs can be cancelled upon
  pressing Ctrl+C, which is not possible with the generic --cluster
  support.

And here:

  If available, DRMAA is preferable over the generic cluster modes
  because it provides better control and error handling.

So I have a conceptual understanding of the advantages of using --drmaa. However, I don't consider the above a very complete explanation, and I don't know how these flags are implemented in snakemake under the hood, can someone elaborate?

Note, that although this could be considered a more general programming question, snakemake is predominantly used in bioinformatics; I was persuaded that this question would be considered on-topic because of this meta-post and this answer.

I'd always kind of wondered how this worked too, so I took this as an excuse to look into the snakemake code. At the end of the day this becomes a question of (1) how are jobs actually submitted and (2) how is it determined if jobs are done (and then whether they failed)?

For DRMAA, python has a module (appropriately named "drmaa") that wraps around the libdrmaa library that comes with most schedulers. This is a very popular route, for example the Galaxy project uses this for dealing with most clusters (e.g., I use it to connect our internal Galaxy instance to our slurm cluster). The huge benefit here is that DRMAA does magic to allow you to simply submit commands to your cluster without needing to know whether you should execute qsub or srun or something else. Further, it then provides methods to simply query whether a job is running or not and what its exit status was.

Using the --cluster command requires a LOT more magic on the side of snakemake. At the end of the day it creates a shell script that it then submits using the command you provided. Importantly, it includes some secret files in that script that it can then watch for (ever noticed that .snakemake directory where you execute a script? That appears to be among the things it's used for.). These are named {jobid}.jobfinished and {jobid}.jobfailed and one of them will get touched, depending on the exit status of your command/script. Once one of those is there, then snakemake can move on in its DAG (or not, if there's a failure). This is obviously a LOT more to keep track of and it's then not possible for snakemake to cancel running jobs, which is something one can do easily with DRMAA.


BigWig circular plot representation I would like to ask if anyone knows or has experience in representing BigWig signals in a circular plot? I mean like IGV but in circular plot not linear. Or any library that can be used in plotting BigWig signals especially in Javascript?

I appreciate any suggestion,
Thank you.

circosJS should probably work well, based on your requirements: https://github.com/nicgirault/circosJS
I posted a short Observable notebook that shows what calls could be used to bring data from remote bigWig files into a circosJS plot:
https://observablehq.com/@alexpreynolds/rendering-bigwig-file-data-in-a-circular-representation
It is less intended as an explicit how-to in Javascript, than as a starting point for the end user with some of the function calls and data scale issues involved. Hopefully it can be useful to the reader for exploratory purposes.


What is the difference between ENCODE Tier 1, 2 and 3 cell types? The ENCODE Experiment Matrix at UCSC lists the different available cell types under the categories "Tier 1", "Tier 2" and "Tier 3". What is the difference between these classifications?

What, for example, makes GM12878 a Tier 1 cell type and A549 a Tier 2 cell type?

The tiers denote the original priority of sequencing and processing the samples. This is understandable given the number of marks and experiment types that ENCODE tried to sequence. Further details are available on ENCODE's website.

Edit: Just to expand a bit, tier one cells were supposed to get processed first. Then tier 2 came second. Of course tier 2 had a LOT more cells, so it got split into 2 and 2.5, since there were already apparently tier 3 cells.


ambiguous calls in nanopolish I am converting scores into frequencies (from nanopolish results) using Nanopore scripts.
The script filters out lots of rows on the premise that reads where LLR < T * N are ambiguous.
where:
LLR = log-likelihood ratio (log-likelihood methylated / log-likelihood unmethylated)
T   = arbitrary threshold. Above that value the C is methylated.
N   = number of Cs (CpG) in the fragment analyzed.

What? Why? Where?

Assuming that it is the postprocessing script calculate_methylation_frequency.py that you are using, user Wouter de Coster's comment is correct that the methylation event is filtered out at the frequency estimation stage if it has a low LLR. It will have  a low LLR if it is considered not confident enough to report. This means not only LLR > T * N, but also the statistical threshold LLR > 0 (see below).
Specifically, here is the code in the script:

def update_call_stats(key, num_called_cpg_sites, is_methylated, sequence):
    if key not in sites:
        sites[key] = SiteStats(num_called_cpg_sites, sequence)

    sites[key].num_reads += 1
    sites[key].called_sites += num_called_cpg_sites
    if is_methylated > 0:
        sites[key].called_sites_methylated += num_called_cpg_sites

That is the function that decides whether to add a site to output, which depends on is_methylated, which is set here:
    # [skip several lines of code]
        num_sites = int(record['num_motifs'])
        llr = float(record['log_lik_ratio'])

        # Skip ambiguous call
        if abs(llr) < args.call_threshold * num_sites:
            # [MP NOTE: SKIP if LLR < T * N]
            continue
        sequence = record['sequence']

        # [MP NOTE: SKIP IF LLR <= 0]
        is_methylated = llr > 0

        # [skip several lines of code]

        update_call_stats(key, num_sites, is_methylated, sequence)

It's a little messy to read because the sites dictionary is a global outside of update_call_stats().
So overall, it seems like if you want sites included that are below this threshold, you will have to either:

use an earlier stage of processing (output of nanopolish call-methylation I believe according to the docs) or
alter the script to change the heuristic thresholds.


Convert Reactome Protein IDs to Pathway IDs? I have a list of about 50k 'Protein IDs' from Reactome. Is there a simple way to get all the corresponding 'Pathway IDs' for each protein? What is the best service to use? (I'm guessing I can use the Reactome API, but I don't necessarily want to hit that 50k times...).

Protein IDs from reactome look like this:

R-HSA-49155
R-HSA-199420

The corresponding Reactome 'Pathway IDs' for those Protein IDs would be:

R-HSA-49155  R-HSA-110331
R-HSA-49155  R-HSA-110330
R-HSA-49155  R-HSA-110357
R-HSA-199420 R-HSA-1660499
R-HSA-199420 R-HSA-202424
R-HSA-199420 R-HSA-199418

If you don't mind hitting it 50k times and are OK with python3...

from urllib import request
import json

def getPathways(proteinID):
    baseURL = 'http://reactome.org/ContentService/data/query'
    PathwayIDs = set()
    try:
        response = request.urlopen('{}/{}'.format(baseURL, proteinID)).read().decode()
        data = json.loads(response)
        if 'consumedByEvent' in data:
            for event in data['consumedByEvent']:
                PathwayIDs.add(event['stId'])
        if 'producedByEvent' in data:
            for event in data['producedByEvent']:
                PathwayIDs.add(event['stId'])
    except:
        pass
    return PathwayIDs

Usage would then be something like:

l = ['R-HSA-49155', 'R-HSA-199420', '']
for rid in l:
    ids = getPathways(rid)
    for _ in ids:
        print("{}\t{}".format(rid, _))

Which would produce:

R-HSA-49155 R-HSA-110239
R-HSA-49155 R-HSA-110240
R-HSA-49155 R-HSA-110238
R-HSA-49155 R-HSA-110356
R-HSA-199420    R-HSA-8948800
R-HSA-199420    R-HSA-6807106
R-HSA-199420    R-HSA-6807206
R-HSA-199420    R-HSA-6807126
R-HSA-199420    R-HSA-8847968
R-HSA-199420    R-HSA-8850997
R-HSA-199420    R-HSA-2321904
R-HSA-199420    R-HSA-8948775
R-HSA-199420    R-HSA-8944497
R-HSA-199420    R-HSA-6807134
R-HSA-199420    R-HSA-8850945
R-HSA-199420    R-HSA-8873946

Note that this will silently ignore invalid or missing IDs such as '' (that's the try and except above. Note also that these are different pathway IDs than what you provided in your example. The main reason is that the protein IDs you showed are not always involved in the pathways IDs you showed (in my example, they always are).


How are MACS2's narrow peak and broad peak algorithms different? The peak calling tool MACS2 can call peaks in either narrow peak mode (for focused signals like transcription factor ChIPseq) or broad peak mode (for more defuse signals, like certain histone modifications).

The algorithm for narrow peak calling is well described in the MACS publication. But I don't find much documentation for how peak calling is different in broad peak mode. The manual only contains the following:

  --broad

  When this flag is on, MACS will try to composite broad regions in
  BED12 ( a gene-model-like format ) by putting nearby highly enriched
  regions into a broad region with loose cutoff. The broad region is
  controlled by another cutoff through --broad-cutoff. The maximum
  length of broad region length is 4 times of d from MACS

But this doesn't really describe exactly how this is performed.

So what is the algorithm MACS uses for calling broad peaks?

The key function is call_broadpeaks:

The description attached to the function says:

  This function try to find enriched regions within which, scores are
  continuously higher than a given cutoff for level 1, and link them
  using the gap above level 2 cutoff with a maximum length of
  lvl2_max_gap.

  scoring_function_s: symbols of functions to calculate
  score. 'p' for pscore, 'q' for qscore, 'f' for fold change, 's' for
  subtraction. for example: ['p', 'q']

  lvl1_cutoff_s:  list of cutoffs
  at highly enriched regions, corresponding to scoring functions.

  lvl2_cutoff_s:  list of cutoffs at less enriched regions,
  corresponding to scoring functions.

  min_length :minimum peak length,
  default 200.

  lvl1_max_gap   :  maximum gap to merge nearby enriched
  peaks, default 50.

  lvl2_max_gap   :  maximum length of linkage
  regions, default 400.

  Return both general PeakIO object for highly
  enriched regions and gapped broad regions in BroadPeakIO.

To give some basic explanation, the algorithm (briefly) appears to be as follows:

Two separate levels of peaks are called, level 1 (a higher pval, ie more significant) and level 2 (a lower pval). Level 1 is controlled by -p and level 2 is controlled by --broad-cutoff. When each peakset is called, they are immediately linked by the max gap parameter for each set.
Then, assuming that all level 1 peaks should be inside level 2 peaks (this is an explicit assumption by MACS2), the algorithm groups level 1 peaks inside level 2 peaks to output a broad peak.

...

This has a few implications:

The broad peak calls really come from the level 2 peaks alone (+ linking). The level 1 peak calls allow you to distinguish sub peaks (so that you can have gapped peaks).
Aside from the linking, the broad peak calls would be the same as narrow peak calls, if you called both with the same pval threshold (for example, if you set --broad-cutoff 0.1 in broad peak mode, and the -p 0.1 for narrow peak mode)


How to filter a SAM file by a bed file? I have a large BAM file and I want to filter it by just what's defined in a bed file, and write it to a new file.

Is it possible to do this with samtools, and if so, how? Is there a better tool I could use to do this?

Thanks to the helpful comments I figured it out:

samtools view -b -h -L bedfile.bed originalbam.bam > newbam.bam

-b = Output as bed file
-h = include header
-L = only output alignments overlapping the input bed file


Looking for genes that are expressed the same in all the samples I have 5 samples (from a RNA seq experiment) from patients but I have no control, so since we can't look for the differentially expressed genes we want to look for the genes that have similar expression in all the sample.
We want to do that to see if some genes that can be involved are among those expressed in all the samples.
I have never done that, I do not know how to do it and I cant find any info about it, so any ideas are welcome.
Thanks!

I cannot think of a way without performing differential expression analysis. Since you don't have controls, what you could do is assigning some of the samples as controls.
Let's say you have samples A, B, C, D, E and F. You can perform differential expression analysis A, B, C (i.e. treatment) versus D, E (control). If you perform this for all possible combinations (that would be 5 choose 2 which is equal to 5 choose 3 if you decide to go for 3 vs 2 splits) and filter for genes which do not show up in any DE list you end up with the genes that are expressed similarly in al of these samples.


Short reads vs long reads with regards to break points My professor in class said the following statements.

Short reads can map break points
Long reads can disambiguate repeat related issues

I'm very new to genomics, so I don't quite understand the first statement. What is meant by short reads can map break points. What are break points to begin with?

A "break point" typically denotes one end of a structural variation, such as an inversion or deletion. You can simply think of it as the position where the normal chromosome was broken to allow the variation. One can find these with short (Illumina) reads relatively easily since with paired end sequencing you'll have pairs spanning these breakpoints with odd relative orientations or insert sizes or such. Note that this is not a particular strength of short reads, since long reads are probably even better at this (they are not as good for finding SNPs at the moment though).


insert size pre and post trimming I have a problem here with my rna seq data:
Sequencing details:
rRNA was removed, followed by cDNA preparation and generation of stranded libraries using the TruSeq Stranded Total RNA Sample Prep Kit. Sequencing performed on the HiSeq2500 platform (Illumina) to generate 2 × 125 bp paired-end reads
Alignment and preprocessing
reads were aligned using tophat2(std aligner in pipeline at core (--library-type fr-firststrand), unfortunately the original unaligned files were purged and i only have access to this aligned file which i converted to fastq using following steps

samtools sort -n sample.bam -o sample_sorted.bam
bedtools bamtofastq -i sample_sorted.bam -fq sample_1.fq -fq2 sample_2.fq (get a lot of mate skipping errors here)
check for these reads for adapters and trim them and again use hisat2 (with --rna-strandness RF)

trimmed reads
Reads were trimmed by passing parameters --adapters adapters.file --adapter-trim-end RIGHT --length-dist --threads 12 --adapter-min-overlap 7 --max-uncalled 250 --min-read-length 25 -- to FLEXBAR version 2.4
adapters for all samples using multiqc

fastqc sequence length pre trimming

fastqc sequence length post trimming

rmats error

Incorrect readLength. sample.bam has a read length of 114, while readLength param is 125

picard error on aligned bam from hisat2)

ProcessExecutor   2 "Not creating insert size PDF as there are duplicated header names: All_Reads"

snapshot of picard insertsize metrics file

Questions

Is my approach flawed at some step?
why do i get different read lengths error for each sample when i process it through rmats when i specify len as 125?
why is no histogram being produced when using picard insertsize metrics?
is normal to see the graph shift to the right after trimming?

insertsize for all samples (data from picard collectinsertsize complied usingmultiqc)

It's not an unreasonable approach. rMATs is rather picky about its input, but it seems you've noticed that.
rMATs can't handle trimmed reads. It also can't handle soft-clipped reads. You might as well not trim. You'll get a lower alignment rate, but your only other choice would be to either exclude the trimmed reads or trim everything down to the same length.
You have a mixture of RF and FR alignments and the plotting script can't handle that.
Assuming you mean the read length distribution, it's actually shifted to the left, not the right. You had a single length before and a long left tail after trimming.


What is the state of the art for GWAS in terms of the statistical algorithm for either Case/control and Quantitative traits? This question was also asked on Biostars

I'm trying to understand what is the best algorithm for GWAS nowadays.
I know we have many tools available like Plink and Hail, but currently, what is the best algorithm if I won't use any them? Let's say, write down a script in R or Python from scratch. Which statistical algorithm should I use?
Is it Linear mixed models (LMMs)? I'm confused as we can have binary phenotypes (case/control) or quantitative phenotypes. LMM seems to address quantitative ones, but can it be used for case/control as well?
Actually, what is the state of the art for both/each of them? Pair-reviewed papers as references will be appreciated.
Thanks!

I've got the answer to my question at https://www.biostars.org/p/410466/ by chrchang522

"The main regression executed by Plink was introduced by EIGENSTRAT in ~2006; see https://www.nature.com/articles/ng1847 . This is actually straightforward to write in R/Python from scratch; the harder part is optimizing the implementation for large datasets.
The Firth regression added to Plink 2.0 to improve handling of rare variants and imbalanced binary phenotypes was motivated by https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4049324/ .
Mixed linear models provide better statistical power when you have lots of close relatives in your dataset, but are much trickier to solve; actually, this is still a significant research area. Two tools covering parts of the current state-of-the-art are SAIGE (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6119127/ ; handles imbalanced binary phenotypes, but relatively slow) and fastGWA (https://www.nature.com/articles/s41588-019-0530-8 ; great speed, but doesn't support dosage data yet and uses a misspecified model for binary phenotypes)."


Why is a PacBio read length larger than the aligned reference region? I recently had some Iso-Seq sequencing done on my organism catfish on the new Sequel platform and got weird alignments for a size selected 4 Kilobase and up fraction after running the isoseq3 pipeline.  When I investigated the exon intron junctions using data from (Pac Bio Iso-Seq, Illumina RNA-seq) against our private reference genome, I get PacBio Iso-Seq reads where the sequence length is greater than the seq length on the reference for my bam alignments:

Seq length: 6751

Seq length on ref: 5267

When blasting the sequence in NCBI it shows 5% similarity to the publicly-available genome.  I'll attach a JBrowse image with the same RNA sequenced with Illumina paired RNA-seq.  I also included a image and circled the pacbio contig and Illumina is above.

Question 1:

Is this a problem with JBrowse?

Question 2:

Does this mean that some DNA was sequenced with the catfish cDNA?

Question 3:

Why is this aligned sequence's length longer than the length of the reference genome region to which it is aligned?

Thanks

Q1
Probably not - this is probably just displaying the contents of the input (I presume a bam file) exactly as they are encoded.
Q2
The splicing looks a little weird, but the exon-mapping pattern in the 5' part of the gene makes it seem like RNA rather than DNA. I think that the large intron is what makes you think that it is DNA, correct?

Explanations for weird splicing

The weird read is a rare splice variant but a true mRNA.
The weird read is not completely matured and has an intron that has not yet been spliced. The IsoSeq library prep protocol picked this up somehow and it ended up being sequenced.

Q3

Explanations for why Seq length: 6751 and Seq length on ref: 5267

There is probably a large amount of clipping on the 5' or 3' end of the read. This could be polyA or some other transcribed-but-not-yet-completely-matured seq. Or, the extra sequence could be sequencing error, PacBio adapters, or a chimeric molecule.

things to try
How to check what is going on:
I'm guessing that you had to upload a bam file to JBrowse to make the IsoSeq track. Use the sequence name that you got from Jbrowse and run the following command to extract the alignment info from the bam file:
samtools view /path/to/your/isoseq/file.bam | grep 'weird_read_id'

Look at the cigar string to see what it says about 5' or 3' clipping and corroborate with the sequence data.


IonTorrent variant caller failing, lack of ZM: tags in BAM files I have a lot of IonTorrent sequencing data in BAM format that I want to call variants on, using the IonTorrent variant caller. However, for some of the BAM files (for about 12 samples out of ~150), the variant caller stops due to what it states is lack of ZM: tags in a few individual reads.

The thing is, most of the reads in the BAM file do not have ZM: tags (when viewed through samtools view), but these do not apparently cause problems.

What settings do I need to change to get these BAM files to process, or do I just filter out the problematic reads and analyse the remainder?

It turns out, after speaking to tech support, that these are old files from a previous version of the software so won't run through variant caller.

The files themselves need reanalysing from the initial FASTQ -> BAM with the newest version.


Which taxonomy data type to use for abundance graphs and statistical analyses? At our institution we have had some debate about this topic. So let me give some background:

Qiime2 produces a OTU table which shows the number of sequences that matched to a specific OTU. It also produces an Taxonomy table which links all the OTUs to their corresponding taxonomic identity.

So what we were discussing was whether to use the number of sequence hits per OTU for alpha and beta diversity, as well as abundance graphs. Or should we use the Taxonomy table and count the number of occurrences for each taxa.

An example of this would be something like the following:

Here are two OTUs:

D_0__Bacteria;D_1__Cyanobacteria;D_2__Oxyphotobacteria;D_3__Chloroplast;D_4__Virgulinella fragilis;D_5__Virgulinella fragilis;D_6__Virgulinella fragilis and D_0__Bacteria;D_1__Cyanobacteria;D_2__Oxyphotobacteria;D_3__Nostocales;D_4__Oscillatoriaceae;D_5__Oscillatoria PCC-10802;Ambiguous_taxa

If we used the first approach then there would be, for example, 3000 sequence matches for the first OTU and 1500 for the second OTU. However, they represent two different species from the same Class. So abundance wise, the first OTU could be represented twice as much as the second. With the second approach they would share the same abundance percentage.

What is the correct approach for these types of analyses according to the bioinformatics community?

The fundamental answer to this question is the role of allele frequency, i.e. population genetics, within taxonomy. The answer of course is that population genetics and population dynamics should be a integral part of taxonomic designation PROVIDING it is not subject to unreported sample bias. Thus it has to be option 1.

However, your definitions I'm not so sure about because the counts per OTU using your literal definition means every SNP within a species. I don't believe that is what you mean, but that is what you have said. Obviously at a taxonomic focus of a single SNP this isn't  meaningful because there would be as many OTUs as there are sequences and the species definition would be meaningless.

Both approaches suffer from lacking a clear definition to make a species assignment and this is controversial. If you are using a tree to assign species membership, e.g. 16S at 5% identity, then that would generally be acceptable, but some would question it.

You can of course use deep learning algorithms to do this, but this is only as accurate your training and at present (check Bioinformatics journal) the accuracy is at best 90%, which is very low by comparison to phylogeny.

Thus I think it is a non-question, because allelic diversity is extremely important in vast swathes of biological analysis.  Taxonomy is not a disapline in its own right, but a foundational step biologists use to understand the ecotope, ecoregion etc... The only  justification for pursing approach 2 is if there is persistent concern in the bias of the sampling process, to the point where the alleleic (species) frequency becomes meaningless.


Is it okay to use deeptools bamCompare (SES normalization) for comparisons across different ATAC-Seq datasets? We are trying to use deeptools for analysis of ATAC Seq datasets. We have datasets with different sequencing depths and are wondering if bamCompare's SES based normalization is appropriate for comparing two ATACSeq datasets as opposed to ChipSeq datasets. If such normalization is okay to perform, does it make sense to compare difference in values across comparisons? Let's say I do bamCompare on A and B using SES, and then I do bamCompare (using SES) on A and C. Will it then make sense to compare the two comparisons (A-B & A-C)?

If such normalization isn't okay, what are the best tools to normalize ATAC-Seq datasets with different sequencing depths so that one can perform comparisons like above?
Thanks for any suggestions or comments.

I would be very hesitant to use SES normalization in a case like this. We recommend that people run plotFingerprint first and see if there is good separation between the samples and only then using SES. I would not expect there to be drastic separation between the samples for ATAC-seq data, so I don't think it will work well. Our standard ATAC-seq normalization is RPGC (aka, "1x normalization"). That should work well for your use-case as well.

As an aside, if you really want to make statistically-backed comparisons between groups I can recommend the CSAW package.


CRISPR/Cas9 screen analysis with Mageck: paired-end sequencing I want to analyse data from a CRISPR/Cas9 screen (control vs. treatment) and I'm using Mageck (https://sourceforge.net/projects/mageck/). The problem is that I'm working with paired-end sequencing data, do I have to analyse R1 and R2 fastq files separately or is possible to do it at same time with Mageck?

The other answer is incorrect. The original poster is asking about paired-end sequencing data, not paired treatment/control experimental design.
Here is how you analyze paired-end sequencing data (this was an answer I gave on Biostars last year):
https://www.biostars.org/p/303656/


How to define an outgroup to build a robust amino acid tree I building a robust amino acid phylogeny with adequate robustness support (see previous post here).
This is a brief summary of what I have done:

I performed a BLASTp analysis using a protein sequence from species A as query against the full NCBI protein database,
downloaded all the hits in fasta format and I added the sequences I identified with the query sequence in my target genome, which is closely related to species A.
I then aligned the sequences using the Compute - Alignment - MUSCLE option in Seqotron and I saved my alignment file in Phylip format.
I tried to run this command line in RAxML to build a robust maximum likelihood tree:

raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s ~/Desktop/alignment_file.phy
However, the program showed the following error:
Error, the outgroup name "outgroup-p" you specified can not be found in the alignment, exiting ....
Question how can I build an outgroup for my phylogenetic tree, and how do I include it in my alignment?
I appreciate your help and suggestions!

I think you have a bug

raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /Users)/username/Desktop/alignment_file.phy

Try the above, where username is your username. I'll check my codes later. The most important bit is the space between the -p and word outgroup. You can omit the -p 10000 it is only needed  for very geeky maximum likelihood and sets the random number stream.


In SBOL, is there a way to distinguish biological replicates? I have a library of biological implementations (i.e. plasmids) of genetic circuits and their experimental data associated. I want to implement an in silico version with SBOL, but to be consistent with what I have I'm my freezer, I need to identify not only the genetic circuit but also the biological replicate (i.e. the different colonies picked), as they have their own associated experimental data.

Support for distinguishing replicates was added in SBOL 2.2, along with other support for linking data throughout the experiment cycle.  I will explain this in terms of SBOL 3, which is simpler, then give the SBOL 2 translation at the end.
For SBOL 3, all designs, including those for a plasmid, a strain, etc, are specified in terms of the Component class. This is just a "paper" design, however, and may or may not actually exist. When you have a real physical aliquot that you want to talk about, this is represented by an Implementation object: this is your biological replicate.
The Implementation doesn't really have any content of its own; its nature is defined by linking back to the intended design with a prov:wasDerivedFrom link. If you've verified that what you got is what you intended, you can supplement that with a built link as well.
One additional note, however: when you're representing a colony, you're probably not actually representing just a plasmid, but rather are representing a particular strain that has been transformed by a plasmid. The actual data you have may then be one level further removed, since the data has been gathered from a sample under particular growth conditions. SBOL can represent that as well, by describing the sample as a mixture of particular media, inducers, etc., with a transformed strain, which is itself a cell type containing your plasmid.
Here is an example diagram showing just this sort of relationship, of data associated with three replicates of a sample.

You can find this example in the SBOL3 tutorial from COMBINE 2020.
Finally, in SBOL 2, this works exactly the same way, except that instead of a Component your implementation points to a ComponentDefinition or a ModuleDefinition, depending on the nature of the design you're representing. I won't go into that distinction, because it was a bad idea and there's a reason we got rid of it for SBOL 3.


Within and between sample count normalization recently I came across a situation where RNASeq sample quality was expresses trough DESeq sizeFactor. So authors reported quality of their samples with respect to some publicly available datasets as good/bad if the scaling factor of their samples was equal/higher/lower than 1. What is more, they computed it using FPKM's from Cufflinks and executed geometric mean calculation over genes (as the one done in DESeq after log transforming the raw counts). Underlying logic being that with FPKM the samples were normalized "within" and with the additional DESeq, normalization "between", thus ready ready for any downstream analysis (differential expression, AS, etc.)
I haven't seen something like that before nor the use of DESeq scaling factor for estimating the quality of a sample (all samples are from the same species/tissue/condition/etc. )
Can anyone comment on why is this ok/not ok because I cannot wrap my head around the logic (then again I am just trying to replicate the study, no value contribution on my end.)
Thank you

The approach you are describing seems very strange.
Crucially, the Vignette for DESeq2 states that the model only works correctly with unnormalized counts as input:

It is important to provide count matrices as input for DESeq2’s
statistical model (Love, Huber, and Anders 2014) to hold, as only the
count values allow assessing the measurement precision correctly. The
DESeq2 model internally corrects for library size, so transformed or
normalized values such as counts scaled by library size should not be
used as input.

Normalized counts like FPKM are not supported and should be converted to counts using the appropriate import funcion.
I don't really understand the use of the normalization factors for quality control because the number of reads in a sample does not correlate well with the quality. Maybe you could link to the paper in case there is something hidden in the experimental design that would explain this.
If you are interested in quality control a useful metric in my experience is the 3' sequencing bias that can be calculated using, for example, RNA-SeQC.


How can I systematically detect unknown barcode/adapter sequences within a set of samples? I have often downloaded datasets from the SRA where the authors failed to mention which adapters were trimmed during the processing.

Local alignments tend to overcome this obstacle, but it feels a bit barbaric.

fastQC works occasionally to pick them up, but sometimes fails to find the actual adapter sequences.

Usually, I ended up looking up the kits they used and trying to grep for all the possible barcodes.

Is there a more robust/efficient way to do this?

You mention that FastQC "fails to find the actual adapter sequences" - I guess you mean in the Adapter Sequence Contamination plot. However, the kmer and Sequence Content Plots are often useful even when the former fails. I've used these in the past - you can sometimes just read off the adapter sequence from the start of the Sequence Content Plot (or at least see how many bases to trim).


Cufflinks Error: sort order of reads in BAMs must be the same I am running Cufflinks for transcriptome assembly using the .bam file generated by Hisat2. I tried both bam and sorted bam files

cufflinks --no-update-check -o ./ -p 15 accepted_hits.bam

it gives the following error:

Error: sort order of reads in BAMs must be the same

UPDATE: It is working fine now. It seems there were some issues with the reference file. The BAM files generated from Tophat2 also gave the same error. I downloaded the updated reference file and it worked great.

Just add --dta-cufflinks to your Hisat2 command, so that the output alignment file provides the attributes needed by Cufflinks (XS flags). This should do the job.

From manual:

  Report alignments tailored specifically for Cufflinks. In addition to
  what HISAT2 does with the above option (--dta), With this option,
  HISAT2 looks for novel splice sites with three signals (GT/AG, GC/AG,
  AT/AC), but all user-provided splice sites are used irrespective of
  their signals. HISAT2 produces an optional field, XS:A:[+-], for every
  spliced alignment.


How do I differentiate outliers from in-group variations in a DESeq2 PCA plot of 9 samples distributed into 3 conditions? I have a PCA plot from DESeq2's plotPCA(vsd, intgroup=c("conditions")) function. I have 9 samples distributed in to 3 groups of 3 biological replicates each. My reason for trying to visualise all on the same PCA plot is to see that samples from the same group (condition) cluster together and to see a separation between the groups.

The PCA plot I have is does not appear as I expect:

1-Can I just remove a sample which I deem is the outlier and perform a reclustering of the rest on the PCA?

2-What parameter can I use to declare that a point is an outlier rather than an in-group Biological variation?

3-I also wish to find out if it is allowed to plot more than two groups of sample at a time on a PCA plot.

Thank you in advance for your kind hints.

Since you have so few points here, it is hard to tell which are outliers by just looking at the PCA plot.
I suggest you use correlation heatmap for your gene expression and check whether there are obvious outliers.

As for using more than one groups when plotting PCA, I think you can just try passing a vector with many groups to the argument intgroup.


Can structural - equation modeling be used in Bioinformatics? I am a bioinformatics student and I was learning the basics of structural equation modeling (SEM).

I was wondering, are there any applications of SEM on  bioinformatics? Are there any strong sides for use of structural equation models using data on bioinformatics?

I did a quick search and found almost nothing. While it seemed like a learning algorithm, I can't say I figured it out completely.

I have came across them in gene regulatory network analysis before. see here and here. However I haven't seen them in other contexts, but that's not to say they don't exist.


Is there a aminoacid mutation that forms similar protein? Is there an amino-acid shift that will 99% of the time end up resulting the same biochemical function/structure and implications? Like Alanine replaced by Leucine will code in 99% of the cases(at any position that first has Alanine) towards proteins that will function the same without big percentual changes in how good it binds or catalyses.

I do have a big amount of mutations. I'd like to filter them out on chance of biological/phenotypical impact, and since some amino acids look like good substitutes by eye (if I look at the molecular structure)

A mutation is most often likely to be either neutral (no effect) or destabiling (protein misfolds). Changes in activity require changes at the active site or the entranceway to it —which are often only a handful of residues, so are very uncommon.

So the answer to your question is most mutations keep the activity, although a large fraction break it. While to get a new activity the process is more fluid and is not boolean.

NB. Not all proteins are enzymes, but for simplicity I will concentrate on enzymes.

Neutrality and core vs. surface

The outer residues of a protein, the surface residues, are very tolerant to mutations. Unless they are part of an binding surface (interface) or the change results in a charge change that is. The core residues on the other hand are more sensitive to mutations.

Destabilisation

One can do a scan in silico or in vitro or in vivo to determine the mutational landscape of a protein, which can take various forms.

This being Bioinformatics stack overflow I'll give an in silico example. The following for example is a heatmap I did some time back with Rosetta pmut_scan, which is a classic tool for this and shows in the y axis the residue it was mutated to and on the x axis the position. A larger change in Gibbs free energy of folding) relative to wild type (∆∆G) means it is destabiling, but anything less than 5 kcal/mol is neutral. As you can see many are neutral, but some half are destabiling —with a caveat for proline, which require a backbone change.

In the in vitro setting one speaks of melting temperature (DSC assays) or turnover number over Michaelis constant (enzyme specific assay). In the in vivo setting one speaks of a fitness cost of a mutation, i.e. the protein is less active so the organism grows less.

Drug resistance

Drugs inhibit enzymes by binding the active site. Therefore resistance is often costly. There are some nice papers that look into the fitness cost and mutational landscape of HIV protease or β-lactamase. In this situation a strong epistasis is seen. That is, a first mutation that confers resistance to the drug greatly destabilises the protein, while a second compensates for the first mutation.

Change in function

An enzyme can catalyse a reaction well on a given substrate, because it was evolved this way (technically the encoding gene evolved, but shh). This is the physiological activity. However, enzymes can also catalyse side activities on similar substrates with a lower turnover and higher Michaelis constant. This is especially true for substrates it was not selected against. This is called "enzyme promiscuity". The tolerance of the physiological activity to mutation is called "robustness", while the ability to increase promiscuous activities is "plasticity". Different enzymes have a different balance of the two.

So a mutation may change the strength of the physiological activity, but may result in a new promiscuous activity that can be used as a stepping stone for selection to increase this activity. Therefore, as you can see it is not a binary switch, except for a few very specific case.

Best amino acid for changes?

Proline is a very strange amino acid as it is a secondary amino acid. It lacks a hydrogen on the backbone nitrogen so it cannot form the required hydrogen bond for an α-helix. And can adopt a cis-conformation with minimal penalty (whereas it is about 5 kcal/mol for other chiral amino acids). So changes to this amino acid are often very destabilising.

Protein cores tend to be hydrophobic so changes to the hydrophobic residues are often deleterious, especially when going from a small residue to a larger one.

The remaining residues are frequently found on the surface so may be neutral. However, some are found in the active site and may be able to shift the parameters of the ligand space for that enzyme. A change in charge may result in an increased specificity for a different set of substrates for example. A size change from a larger amino acid to a smaller one may result in a broader range of promiscuous activities —mutations to alanine for example.

So ironically, polar and charged residues as they are frequently found on surfaces are the residues which are most often neutral, but in a small amount of cases they may change activity.


Pymol: color spheres with different colors Hello I am a Pymol user
I want to color all alpha carbons (ca) and beta carbons (cb) from a certain ressidue in different color. How can I do this?
I selected the spheres with:
show spheres, C/1-9/ca
show spheres, C/1-9/cb

and if I use
set sphere_color, theColor

I color all the spheres what I don't want.
according to: color    gold,    resi    200-‐221
I tried
gold,    spheres, C/1-9/ca
but this didn't work

Spheres
So how PyMOL stores how it shows atoms is that each atom has a colour and a represention, which is a binary number telling what to show. Therefore the colour for spheres is the same as the colour for sticks. The only way round this is to copy the object (create copycat, original and you set copycat differently)
To set the atom you do
select iliketheschaps, resi 200‐221 and name CA+CB
color tv_green, iliketheschaps
show spheres, iliketheschaps

The object iliketheschaps is a selection. Do make sure the range is a hyphen-minus (character next to 0) as opposed to an en-dash (–) which Word subs and on a mac you get with cmd+-
non-PyMOL colours
Now you asked for "gold". This is not a colour in the menu —famous for the colour limon. Gold is a CSS colour name, but it won't work.
You will have to colour by hex code. Say you want "gold" you need to google the colour name and hex. #FFD700.
color 0xFFD700, iliketheschaps

The Fashion Police have been alerted

This paragraph is in jest.

However, as you may know bioinformaticians are very fashion savvy for colours. Mustard was the in-colour from last year, so you really ought to use mustard (0xffdb58). This year is all about nectarine (0xffa168) (which is not salmon —this ain't 2015!). Coral (0xff7f50) and teal (0x008080) were from previous years and are still are going strong!


Find paralogs in a draft genome We generated a (diploid, chordata, highly heterozgous) genome using PacBio and we wanted to see whether it contains lineage-specific duplications (paralogs, basically). The genome is not in Ensembl yet.

The only data we have at the moment are:

genome
transcript annotation
RNAseq

We found some methods from papers:

use Blast
detect segmental duplications in complete genomes with SDDdetector
detecting putative recent segmental duplications or diploid homozygous organisms based on NGS data DuplicationDetector
(I just came up with this, but should work) map back the reads to the assembly & analyze the read depth to detect duplicated portions

I'll gladly take advice.

This is really tough to do with highly heterozygous animals. What are your genome assembly stats? Specifically, what is your number of contigs, scaffolds, assembly size, and the N50?

If you have proximity ligation data it will be easier to determine if potential paralogs are truly from different regions of the genome or are just from homologous locations on sister chromatids that ended up in uncollapsed in your final assembly.

If I were you and had a good genome, I would start by annotating gene models using RNAseq data, then blasting each gene in the model against all genes in the model. This would give you potential paralogs to start looking at more closely.


Seurat with normalized count matrix? I know that in Seurat we have the function CreateSeuratObject from which the analysis starts, but it accepts raw count matrix according to the documentation. I have only the already normalized count matrix, so is there a way to work with Seurat using normalized data? I would try maybe to feed in the normalized data instead of raw to CreateSeuratObject, but I do not know what Seurat is doing behind the scenes, whether it is a good idea. Besides, Seurat provides by default only one log - normalization method, but I may want to normalize the data by myself with various methods and only then start the analysis with Seurat - that is the other reason why I want to find a way to start from normalized data.

This was addressed by the Seurat developers here:

  if you have TPM counts, I suggest you don't use
  Seurat::NormalizeData(), since TPM counts are already normalized for
  sequencing depth and transcript/gene length. Note that
  Seurat::NormalizeData() normalizes the data for sequencing depth, and
  then transforms it to log space. If you have TPM data, you can simply
  manually log transform the gene expression matrix in the object@data
  slot before scaling the data.

And here:

  If you want to run your own normalization function, please place the
  data in log scale before placing it into object@data

And here:

  If you're supplying already normalized values to CreateSeuratObject,
  you should skip the normalization step. You will likely see a warning
  when you run ScaleData (which is checking whether you normalized using
  Seurat). You can ignore this or turn it off by setting check.for.norm
  to FALSE.

This is a lot of quotes, but I found them all to be helpful in different ways.


Using data mining of papers in order to derive genomic connections Are there any examples of researchers data mining articles and papers (e.g., from pubmed or google scholar) in order to derive relations between genes and diseases?

There are many databases that have used publication scraping for oncogenic gene fusions. There are publications for the individual methods they used for their scraping and aggregation.

COSMIC - http://cancer.sanger.ac.uk/cosmic
TICdb - http://www.unav.es/genetica/TICdb/
ChimerDB - http://ercsb.ewha.ac.kr/fusiongene/ -
https://www.ncbi.nlm.nih.gov/pubmed/27899563


What does PCA mean on GWAS I understand what GWAS is and I'm able to perform certain tests with the p-values, etc. But what I am having a hard time wrapping my head around is what PCA on GWAS means.

So let's say I have 100,000 individuals and genotype data for 10 million SNPs. And I'm analyzing this for 10 different diseases. I also have p-value information for each of the 10 diseases for the 10 million SNPs.

At my work, they were throwing around the following phrases. I'm not quite sure what exactly it means.

"Perform PCA on GWAS." Does this mean, I'm performing PCA on the genotype data or p-value data? I've done PCA on classic machine learning problems like "iris." But I'm not able to quite wrap my head around what PCA information tells me in terms of this genetic data.
"PCA on individual diseases vs. PCA on entire GWAS" What exactly is the difference. I realize this is a broken statement, but they said we have two sets of PCs. One on the individual diseases which are unique, and then PCs on the entire GWAS for all 10 diseases.

Sorry if my questions is not properly asked, I realize it's probably heavily abusing notation.

Added clarifications:

So, my matrix has 100k rows and 10 million columns. Then, PC1, for example, is a vector with 100k values, each value corresponding to each individual. What exactly does PC1 value that's associated with the first individual tell me with regards to the PC1 value that's associated with 2nd individual, and so on. My current understanding if you think of a number line that represents PC1, then each of the PC values (corresponding to each individual) represents where on the number line that individual falls.
How exactly does this eliminate/account for variation in population stratification?

From my memory of what a statistician told me, a PCA aims to determine independent linear combinations of variables (i.e. genotypes) that account for the most variation in the dataset. With 10 million SNPs, the vectors describing the genotypes of all individuals reside somewhere inside a 10-million-dimension hypersphere, and the first principle component determines the best vantage point / angle from outside that sphere that will create the greatest spread of data when data is projected onto a perpendicular / orthogonal axis. The second principle component is the perpendicular / orthogonal axis that will create the second-greatest spread of data.

Aside: There's some wonderful linear algebra (which I've forgotten) that makes it very quick (computationally) to do this angle calculation, exploiting the nature of linear transformations and orthogonality.

The point of the PCA for GWAS is to identify (and preferably exclude) differences that might be due to the structure of populations, rather than the disease of interest. It's fairly common for researchers to do a PCA on the entire set of genotype data, then use the first 5, 10, or 20 principal components as covariates in the association model.

There's an underlying assumption that the majority of genetic variation is due to population structure - if this is not the case, then the disease would have to be extremely omnigenetic, with almost the entirety of the genome contributing to the disease. The approach also won't deal with all the population structure (because there may be non-linear contributions of genotypes to population structure), but it will get rid of the largest problems.

Additional to this, I'd recommend the following things:

Structure Visualisation

Carry out a STRUCTURE analysis or PCA on the combined cases and controls, and display the results of the first two principle components, labelling, colouring, or otherwise indicating both groups. As an additional check, it may be useful to do addititional PCAs for cases and controls separately.

If they appear different, it suggests that the cases and controls were drawn from separate population clusters, which should be commented on in the discussion of the paper. There's no need to panic and find more data, just make the observation.

Subsampling

Carry out association tests on an equal-sized subsample of both datasets. This will mean that area-under-the-curve statistics can be properly interpreted. If you've got unbalanced datasets, use a little bit less than the smaller of the two groups (e.g. for 5k diseased and 95k healthy, use a random subsample of 4.5k from both groups).

Multiple testing / Internal validation

Run the association tests multiple times (with different subsamples), recording the individuals that were subsampled for future reference / reproducibility.

If an association is only present in one of the many different subsamples, it should be ignored as a spurious result. It means that there is some association that is not generalisable to all (or most) of the cases and controls. These variants might be interesting in themselves on an individual level, but not in a general population.

I explain some of this in more detail in the following preprint (which I'm not planning to submit to any journal, due to it being a bit stale and needing a lot of work to fix up):

Bootstrap Distillation: Non-parametric Internal Validation of GWAS Results by Subgroup Resampling


Global or local alignment on sequences for which it is assumed that they have common ancestry As the title says I would like to know which alignment method(global or local) should I use for sequences for which is assumed that they have common ancestry?

The sequences are 21 aa long peptide chains (proteins).

About the background of the sequences, it's only stated that it's ASSUMED that they are from the same ancestry.

The thing that we assume that they are from the same ancestry makes me confused. I know that on similar sequences with the same length usually, we use global.

My initial reaction is that 21 amino acids is too little information to demonstrate these proteins are homologues (which is the question behind your question). A homologous site means it shares a common ancestry and is distinguishable from random chance or convergent evolution. A common ancestor can be demonstrated using phylogeny (trees).

If a regular alignment program cannot align these sequences, such as clustalo or muscle, the assumption they are homologues cannot be easily tested using phylogeny. That would then become complicated at which point you would need to ask whether it is worth pursing the question.

My recommendation is to perform a simple protein-protein blast and assess the e-scores to get some feel about you data. Your question is achievable, but you would solve this in context to additional sequence data. Whether that sequence data is a denser sample or else contiguous sequence data needs further exploration.