text:
Detecting broad peaks in sRNA-seq data
What kind of tool would be appropriate do detect "broad peaks" in small RNA-seq sequencing data?
MACS2 appears to be  developed for ChIP-seq data, but I see that there is a --nomodel option. Would that make this program usable in my case?
I suppose I should also set the --shift and --extsize options, but I'm not sure how:
  --nomodel             Whether or not to build the shifting model. If True,
                        MACS will not build model. by default it means
                        shifting size = 100, try to set extsize to change it.
                        DEFAULT: False
  --shift SHIFT         (NOT the legacy --shiftsize option!) The arbitrary
                        shift in bp. Use discretion while setting it other
                        than default value. When NOMODEL is set, MACS will use
                        this value to move cutting ends (5') towards 5'->3'
                        direction then apply EXTSIZE to extend them to
                        fragments. When this value is negative, ends will be
                        moved toward 3'->5' direction. Recommended to keep it
                        as default 0 for ChIP-Seq datasets, or -1 * half of
                        EXTSIZE together with EXTSIZE option for detecting
                        enriched cutting loci such as certain DNAseI-Seq
                        datasets. Note, you can't set values other than 0 if
                        format is BAMPE for paired-end data. DEFAULT: 0.
  --extsize EXTSIZE     The arbitrary extension size in bp. When nomodel is
                        true, MACS will use this value as fragment size to
                        extend each read towards 3' end, then pile them up.
                        It's exactly twice the number of obsolete SHIFTSIZE.
                        In previous language, each read is moved 5'->3'
                        direction to middle of fragment by 1/2 d, then
                        extended to both direction with 1/2 d. This is
                        equivalent to say each read is extended towards 5'->3'
                        into a d size fragment. DEFAULT: 200. EXTSIZE and
                        SHIFT can be combined when necessary. Check SHIFT
                        option

In my case the tags are the same things as the fragments, so this "extension" thing is probably not appropriate.
Basically, I would like to call peaks based on the plain real coverage of my small reads. Are there other tools that would be more appropriate?

Do you actually need peaks or is your end goal instead to find novel sRNAs? I would think stringTie followed by a considerable amount of filtering would be better.

I mapped sRNA-seq reads on the genome, and I want to automatically delineate genomic regions that produce them. By "broad peaks", I mean such regions. Not sure what you mean by "novel sRNAs".

sRNAs shouldn’t appear as broad peaks unless they’re coming from repeat clusters. See the answer from Sebastian Müller for good tool options.

The sure why any software should be confined to only broad peaks, but there is two packages that should do what I thing you have in mind. Just to clarify, for smallRNAs peak detection is sometimes also referred to as loci or cluster detection:

ShortStack which is developed for alignment, annotation, and quantification of small RNAs. The de-novo cluster finding option as copied from the manual might do the trick:

DE-NOVO CLUSTER FINDING
    Unless options --locus or --locifile are used (see below), ShortStack
    will de-novo identify clusters of small RNA accumulation genome-wide.
    Cluster definition is simple: First, all regions containing at least one
    primary alignment are found where the maximum distance between the ends
    of the alignments is <= option --pad (default: 75). Second, if the
    number of alignments in the cluster is >= option --mincov (default:
    0.5rpm), the cluster is kept. The mincov threshold can also be specified
    in terms of reads per million by using a value such as 3.2rpm (which
    specifies the threshold to be 3.2 reads per million). Using a rpm
    threshold allows the sensitivity of cluster discovery to be normalized
    between libraries of different sizes. Alternatively, reads per million
    mapped (rpmm) can be specified: A --mincov of 1.2rpmm indicates 1.2
    reads per million mapped is the threshold. rpm is a fraction of total
    library size, while rpmm is a fraction of only the aligned & placed
    fraction of the library.

segmentSeq
is an R-package useing a Bayesian approach for small RNA locus detection. It's however really slow and not developed any further by the main author.


text:
Automating sequence download from Wormbase
I'm working on the C. elegans model organism, for which there exists an online biological information database: Wormbase.
I'm able to manually search for a transcript (say "K06C4.12") and this gets me the following URL:
http://www.wormbase.org/species/c_elegans/cds/K06C4.12
Say I have this transcript ID: "F08G2.2". I guess that I can form a valid similar URL, as follows:
http://www.wormbase.org/species/c_elegans/cds/F08G2.2
Indeed, it works as expected.
Now, I would like to go further and be able to download in an automated way what can be obtained by clicking on one of the "download" buttons in the "Sequence" section of the above pages. More specifically, I'm interested in the "spliced + UTR" and "conceptual translation" fasta sequences.
If I use the DOM inspector in Firefox (Shift-Ctrl-C), I can locate things like the following:

    download

        >K06C4.12 spliced + UTR
        ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTCACCAAGCCAAAGGACGGAAAGAAGAG
        ACGTCATGCCCGTAAAGAATCCTACTCCGTGTACATCTACCGTGTGCTCAAGCAGGTTCATCCAGACACTGGAGTCTCTT
        CCAAGGCTATGTCCATCATGAACTCTTTCGTCAACGACGTCTTCGAGCGTATCGCTGCTGAGGCATCTCGTCTTGCTCAC
        TACAACAAGAGATCCACCATCTCTTCTCGCGAAATTCAAACTGCCGTCCGTTTGATTCTTCCAGGAGAGCTTGCCAAGCA
        CGCCGTCTCTGAGGGAACCAAGGCTGTCACCAAGTACACCTCTAGCAAGTAA

K06C4.12conceptual translation.fasta

        >conceptual translation
        MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSSKAMSIMNSFVNDVFERIAAEASRLAH
        YNKRSTISSREIQTAVRLILPGELAKHAVSEGTKAVTKYTSSK*

    view conceptual translation (123 aa)

    >conceptual translation

    MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSSKAMSIMNSFVNDVFERIAAEASRLAH

    YNKRSTISSREIQTAVRLILPGELAKHAVSEGTKAVTKYTSSK*

I don't find these if I simply download the page with wget and search in the source code, so they must be dynamically generated.
How could I proceed to get the sequences?
The documentation about Wormbases's REST API is not very detailed, so I don't understand how to form an URL that would give me the desired sequences.

Did you try to download the whole database and search there? (I am not sure if you need large amounts of data or not ...)

I don't know where to download it, neither how it is structured. I don't need a large amount of data, so I hope querying the site should be OK.

To download it they said it is here. There should be a ftp site, so you could try that too.

Indeed, they even provide a scary shell script (first time I see a "goto" instruction in a shell script) to install the database: ftp://ftp.wormbase.org/pub/wormbase/releases/current-production-release/acedb/INSTALL
I'm not sure I'm ready to delve into this for the moment...

Further inspecting the source code of the downloaded page, I found the following:

    <span class="ui-icon ui-icon-close module-close" wname="sequences"

    Sequences

This suggested me to have a look at http://www.wormbase.org/rest/widget/transcript/K06C4.12/sequences and I obtained a page containing the elements found in the DOM inspector.
This can be parsed, for instance using the bs4 python module:
# To get the page
import requests
# To parse the page content
from bs4 import BeautifulSoup
# To get this into Biopython
from io import StringIO
from Bio import SeqIO

rest_template = "http://www.wormbase.org/rest/widget/transcript/{accession}/sequences"
rest = rest_template.format(accession="K06C4.12")
# This did not work with "html.parser"
soup = BeautifulSoup(requests.get(rest).text, "html5lib")

def generate_fasta(soup):
    name = None
    for div in soup.find_all("div"):
        if "ui-helper-hidden" in div["class"]:
            if div["id"] == "content":
                fasta = SeqIO.read(StringIO(div.text), "fasta")
                description_fields = fasta.description.split()
                if len(description_fields) == 4:
                    fasta.id = "_".join(description_fields[:2])
                    if name is None:
                        name = fasta.name
                elif len(description_fields) == 2:
                    # This assumes that the "conceptual_translation"
                    # is after the "unspliced + UTR" or "spliced + UTR"
                    # Otherwise, one would need to pass
                    # the accession as function parameter
                    fasta.id = f"{name}_translation"
                else:
                    fasta.id = "_".join(description_fields)
                yield (fasta.id, fasta)

records_dict = dict(generate_fasta(soup))
print(records_dict['K06C4.12_spliced'])
print(records_dict['K06C4.12_translation'])

Output:
ID: K06C4.12_spliced
Name: K06C4.12
Description: K06C4.12 spliced + UTR
Number of features: 0
Seq('ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTC...TAA', SingleLetterAlphabet())
ID: K06C4.12_translation
Name: conceptual
Description: conceptual translation
Number of features: 0
Seq('MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSS...SK*', SingleLetterAlphabet())

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
PCA vs tSNE in single cell RNA-seq
What makes tSNE being the preferred dimensional reduction for visualization in single cell RNA-seq over PCA?
I am aware that tSNE works better at showing local structures and fails to capture global structures of the data.
But I think I don't fully get the reason of why is this an advantage? It offers better resolution? Or better separation of the cells as compared to PCA?

They are often stacked together, i.e. first PCA then tSNE. They do quite different things.

tSNE often offers better visual representation (separation) on such complicated data than PCA. As Micheal pointed out, computing a tSNE embedding over 20.000 gene dimensions is computationally unfeasible, so a number of PCs are normally calculated and these are used as input for calculating the tSNE. They are used in tandem.
As for global vs. local, we are much more interested to see similarity of cells to a limited number of neighbours indicating a celltype and grouping these close together. This is more important than the distance between such celltypes. (assuming the separation in your tSNE is driven by a biologically meaningful factor such as celltype    and not some confounder)
Edit:
Since I just made these images for myself anyway I might as well post them. The first 2 PCs show some separation by celltype, tSNE computed over 100 PCs gets very nice separation.

UMAP is another dimensionality reduction algorithm that better respects global structure.

Alse checkout sleepwalk which helps you see where your embedding works well and where it doesn't https://anders-biostat.github.io/sleepwalk/

Woooaaa this a quite a dramatic demonstration. Only one upvote per post ... pity

Woaaw! Super clear now!

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Efficiently iterate over bed and subset pyranges objects
I am using pyranges and I have two datasets that I want to compare in several genome-wide intervals based on a bed file.
For every interval in the bed file, I want to get all overlapping positions from the two datasets and sum two columns. That gives me 4 numbers for a fisher exact test, but that's not the most important part.
My current code looks like below, and just iterates over the itertuples() of the bed file as pandas dataframe, and uses that to slice the pyranges object:

import pandas as pd
import pyranges as pr
from pyranges import PyRanges
from scipy.stats import fisher_exact
import numpy as np

def main():
    colnames = ["Chromosome", "Start", "End", "num_motifs", "calls", "methylated", "freq", "seq"]
    h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0))
    h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0))
    for chr, begin, end in pr.read_bed("chr21.bed.gz").merge().as_df().itertuples(index=False, name=None):
        odr, pv = fisher_exact([meth_counts(h1[chr, begin:end]),
                                meth_counts(h2[chr, begin:end])])

def meth_counts(interval):
    if not interval.empty:
        all, methylated = interval.as_df()[['calls', "methylated"]].sum()
        return [all - methylated, methylated]
    else:
        return [0, 0]

This code works well, but I'm not sure if I use the "full power" of what pyranges could offer me for iterating and subsetting. I uploaded some example data on google drive (250GB).

Note that gr[chr, begin:end] will create a new NCLS to do lookups for every single iteration. This is sloooow. It is better to use functions like overlap, intersect, join etc.

Setup:
import pandas as pd
import pyranges as pr
from pyranges import PyRanges
from scipy.stats import fisher_exact
import numpy as np

# ! zcat dataset1.tsv.gz | head -2
# chromosome    start   end num_motifs_in_group called_sites    called_sites_methylated methylated_frequency    group_sequence
# chr21 5010053 5010053 1   3   0   0.000   CACCACGTCCA

# avoid using more columns than we need, especially seq which contains object
# is heavy to lug around
colnames = ["Chromosome", "Start", "End", "calls", "methylated"]
usecols = [0, 1, 2, 4, 5]

h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))
h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))

regions = pr.read_bed("chr21.bed.gz")
regions = regions[~regions.Chromosome.str.contains("_")].merge()

regions.ID = list(range(len(regions)))

regions
# +--------------+-----------+-----------+-----------+
# | Chromosome   | Start     | End       | ID        |
# | (category)   | (int32)   | (int32)   | (int64)   |
# |--------------+-----------+-----------+-----------|
# | chr21        | 5065660   | 5065810   | 0         |
# | chr21        | 5068140   | 5068290   | 1         |
# | chr21        | 5128400   | 5128550   | 2         |
# | chr21        | 5154780   | 5154930   | 3         |
# | ...          | ...       | ...       | ...       |
# | chr21        | 46637040  | 46637190  | 2377      |
# | chr21        | 46638480  | 46638630  | 2378      |
# | chr21        | 46639820  | 46639970  | 2379      |
# | chr21        | 46692680  | 46692830  | 2380      |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 2,381 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

h1
# +--------------+-----------+-----------+-----------+--------------+
# | Chromosome   | Start     | End       | calls     | methylated   |
# | (category)   | (int32)   | (int32)   | (int64)   | (int64)      |
# |--------------+-----------+-----------+-----------+--------------|
# | chr21        | 5010053   | 5010053   | 3         | 0            |
# | chr21        | 5010215   | 5010215   | 5         | 1            |
# | chr21        | 5010331   | 5010335   | 6         | 0            |
# | chr21        | 5010525   | 5010525   | 2         | 2            |
# | ...          | ...       | ...       | ...       | ...          |
# | chr21        | 46699713  | 46699718  | 274       | 254          |
# | chr21        | 46699739  | 46699739  | 133       | 113          |
# | chr21        | 46699756  | 46699756  | 37        | 32           |
# | chr21        | 46699770  | 46699773  | 156       | 130          |
# +--------------+-----------+-----------+-----------+--------------+
# Unstranded PyRanges object has 331,373 rows and 5 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Get total number of methylated per region:
def count_total_and_methylated(df):
    """For each region in bed, get sum of calls and methylated"""

    grpby = df.groupby("ID")
    total = grpby.calls.sum()
    methylated = grpby.methylated.sum()

    chromosome = df.Chromosome.iloc[0]
    start = grpby.Start.first()
    end = grpby.End.first()
    region_id = grpby.ID.first()

    return pd.DataFrame({"Chromosome": chromosome, "Start": start, "End": end, "ID": region_id,
                         "calls": total, "methylated": methylated})

rh1 = regions.join(h1).drop(like="_b").apply(count_total_and_methylated)
rh2 = regions.join(h2).drop(like="_b").apply(count_total_and_methylated)

rh1.File = 1
rh2.File = 2

def add_missing(df, df2):

    """Since join is inner rh1 and rh2 might not have the exact same intervals. Fix this."""

    missing_in_first = np.setdiff1d(df2.ID, df.ID)
    first_missing = df2[df2.ID.isin(missing_in_first)].copy()

    f = df.File.head().iloc[0]

    first_missing.loc[:, "File"] = f
    first_missing.loc[:, "calls"] = 0
    first_missing.loc[:, "methylated"] = 0

    return pd.concat([df, first_missing])

rh1 = rh1.apply_pair(rh2, add_missing)
rh2 = rh2.apply_pair(rh1, add_missing)

j = rh1.join(rh2, suffix=2).drop(like="File*|ID*|(Start|End)2")
j
# +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
# | Chromosome   | Start     | End       | calls     | methylated   | calls2    | methylated2   |
# | (object)     | (int32)   | (int32)   | (int64)   | (int64)      | (int64)   | (int64)       |
# |--------------+-----------+-----------+-----------+--------------+-----------+---------------|
# | chr21        | 5065660   | 5065810   | 188       | 8            | 70        | 0             |
# | chr21        | 5068140   | 5068290   | 168       | 166          | 74        | 72            |
# | chr21        | 5128400   | 5128550   | 901       | 14           | 143       | 1             |
# | chr21        | 5154780   | 5154930   | 653       | 0            | 593       | 0             |
# | ...          | ...       | ...       | ...       | ...          | ...       | ...           |
# | chr21        | 46637040  | 46637190  | 66        | 63           | 69        | 62            |
# | chr21        | 46638480  | 46638630  | 28        | 0            | 30        | 3             |
# | chr21        | 46639820  | 46639970  | 22        | 0            | 18        | 14            |
# | chr21        | 46692680  | 46692830  | 71        | 25           | 57        | 23            |
# +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
# Unstranded PyRanges object has 1,922 rows and 7 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Compute fisher exact using fisher instead of scipy:
def fe(df):
    from fisher import pvalue_npy

    cols = "calls methylated calls2 methylated2".split()
    df[cols] = df[cols].astype(np.uint)
    _, _, two_sided = pvalue_npy(df.calls.values, df.methylated.values, df.calls2.values, df.methylated2.values)

    t1, t2, m1, m2 = df.calls, df.calls2, df.methylated, df.methylated2

    _or =   (t1 * m2) / (t2 * m1)

    df.insert(df.shape[1], "PValue", two_sided)
    df.insert(df.shape[1], "OR", _or)

    return df

result = j.apply(fe, nb_cpu=1)
result.drop(like="methylated*|calls*")
# +--------------+-----------+-----------+------------------------+--------------------+
# | Chromosome   | Start     | End       | PValue                 | OR                 |
# | (object)     | (int32)   | (int32)   | (float64)              | (float64)          |
# |--------------+-----------+-----------+------------------------+--------------------|
# | chr21        | 5065660   | 5065810   | 0.1152954430992567     | 0.0                |
# | chr21        | 5068140   | 5068290   | 0.9999999999995065     | 0.9846955389124064 |
# | chr21        | 5128400   | 5128550   | 0.7075919513120225     | 0.4500499500499501 |
# | chr21        | 5154780   | 5154930   | 1.0                    | nan                |
# | ...          | ...       | ...       | ...                    | ...                |
# | chr21        | 46638480  | 46638630  | 0.24262295081967644    | inf                |
# | chr21        | 46639820  | 46639970  | 0.00024909594620266284 | inf                |
# | chr21        | 46692680  | 46692830  | 0.7354565511951735     | 1.1459649122807019 |
# | chr21        | 36784100  | 36784250  | 1.0                    | nan                |
# +--------------+-----------+-----------+------------------------+--------------------+
# Unstranded PyRanges object has 1,932 rows and 5 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Thanks a lot for your reply. It seems like a decent solution but I need some more time to test it, as I'm focussing on getting a grant written these days...

Yeah, I'm writing a PR now. But with better code than here :)

I've tested that it gives the same results as your code. And it is super fast. Just wait for my PR plz

Unless you need to use Python and can't use subprocess, here's a quick CLI one-liner which sums signal from sorted BED5 files, over the genomic space where they overlap:
$ bedmap --echo --sum --delim '\t' <(bedops --merge A.bed B.bed ... N.bed) <(bedops --everything A.bed B.bed ... N.bed) > answer.bed

(Requires bash for process substitutions.)

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Different line length in sequence 'chrY'
I just downloaded a reference genome (using wget), and attempted to use it with samtools view. However, I received the following error:
[E::fai_build_core] Different line length in sequence 'chrY'
[E::fai_path] Failed to build index file for reference file 'GRCh38_full_analysis_set_plus_decoy_hla.fa'
[E::fai_build_core] Different line length in sequence 'chrY'
[E::refs_load_fai] Failed to open reference file 'GRCh38_full_analysis_set_plus_decoy_hla.fa'
[E::hts_open_format] Failed to open file "ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR324/ERR3240208/HG00244.final.cram" : No such file or directory
samtools view: failed to open "ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR324/ERR3240208/HG00244.final.cram" for reading: No such file or directory

Just looking at the beginning of chrY in the fasta file does not look suspicious (grep chrY ).

The error message is accurate as the downloaded fasta file got corrupted.
Printing all the contigs, using the command below (the quotes around > are important, without quotes you will erase your fasta!)
grep ">" /results/rr/study/hg38s/study252-P200_analysis/fus-analysis/1000G_NYGC/GRCh38_full_analysis_set_plus_decoy_hla.fa | less -S

returns:
>chr20  AC:CM000682.2  gi:568336004  LN:64444167  rl:Chromosome  M5:b18e6c531b0bd70e949a7fc20859cb01  AS:GRCh38
>chr21  AC:CM000683.2  gi:568336003  LN:46709983  rl:Chromosome  M5:974dc7aec0b755b19f031418fdedf293  AS:GRCh38  hm:multiple
>chr22  AC:CM000684.2  gi:568336002  LN:50818468  rl:Chromosome  M5:ac37ec46683600f808cdd41eac1d55cd  AS:GRCh38  hm:multiple
>chrX  AC:CM000685.2  gi:568336001  LN:156040895  rl:Chromosome  M5:2b3a55ff7f58eb308420c8a9b11cac50  AS:GRCh38
>chrY  AC:CM000686.2  gi:568336000  LN:57227415  rl:Chromosome  M5:ce3e31103314a704255f3cd90369ecce  AS:GRCh38  hm:10001-2781479,56887903-57217415
CATCTTGAATTAATTTTTGTGTAAGGTGTAAGGAAGGGATCC<Error>ConnectionClosedException<Message>Premature end of Content-Length delimited message body (expected: 189,622,683; received: 26,672,030)</Message><ErrorMessa>
>chrM  AC:J01415.2  gi:113200490  LN:16569  rl:Mitochondrion  M5:c68f52674c9fb33aef52dcf399755519  AS:GRCh38  tp:circular
>chr1_KI270706v1_random  AC:KI270706.1  gi:568335410  LN:175055  rg:chr1  rl:unlocalized  M5:62def1a794b3e18192863d187af956e6  AS:GRCh38
>chr1_KI270707v1_random  AC:KI270707.1  gi:568335409  LN:32032  rg:chr1  rl:unlocalized  M5:78135804eb15220565483b7cdd02f3be  AS:GRCh38
>chr1_KI270708v1_random  AC:KI270708.1  gi:568335408  LN:127682  rg:chr1  rl:unlocalized  M5:1e95e047b98ed92148dd84d6c037158c  AS:GRCh38
>chr1_KI270709v1_random  AC:KI270709.1  gi:568335407  LN:66860  rg:chr1  rl:unlocalized  M5:4e2db2933ea96aee8dab54af60ecb37d  AS:GRCh38

Showing a download error (in this case) somewhere in chrY.

Thanks for posting the answer! Please remember to accept it so the question can be marked as answered (I know you can't do it immediately, that's why I'm leaving this as a reminder).


text:
Classification (supervised learning) of expression data on pathway level
I was wondering if there is any way to apply classification algorithms (e.g random forest) on microarray data but not using the genes as predictors/features but the pathway they belong to.
The thing is that the expressions of the genes that belong to the same pathway should somehow, grouped together and transformed into a new value that will describe that pathway for each specific sample.
Let me give an example.
Assuming that we have g1, g2, g3, g4, g5, g6, g7, g8, g9, g10 as genes
and p1, p2, p3 as pathways.
And g1, g2, g3 belong to p1 while g3, g4, g5, g6belong to p2 and g2, g6, g7, g8, g9, g10 belong to p3
One thought is to calculate the mean expression of the above gene groups per sample and then run the classification algorithm on top of these values. But I would like to hear/read your opinion on such an approach and also your suggestions on using a different method to do this.
The final goal would be the extraction of the most significant classification predictors (pathways).

Do you assume that expression of genes are similar in pathways? Because in a real experiment they are not.

No, I don't want to assume something like that but I mentioned the mean in order to group the genes of each pathway. This approach might be totally wrong but that's why I posted here. To get feedback and different opinions.

My point is, that I haven't seen an algorithm that can transform the expression values of different genes into one value that represents a given pathway.

Oh ok. I see...

It looks like you're trying to incorrectly reinvent the activation Z-score produced by IPA or similar tools. I suggest you simply use that instead.

IPA is commercial. Isn't it? Can you provide me with free R-based ones to understand what exactly you suggest?

Have a look at the GSVA package. It allows to convert a matrix with genes x Samples to a pathways x Samples using several methods ssgsea, gage, gsva...
Afterwards you can use that matrix as input for differential expression of pathways or classification algorithms or whatever.

However it depends on the input of the "pathways" you give it. Make sure that the pathways are correct. Also consider that the pathway topology might change with conditions.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Finding new drug-protein targets/interactions using a known one as base
Given a known interaction of a molecule with some proteases (having X-rays images too) is it possible to predict interactions with other proteins either using their 3D structure or not?
If anyone has already tried such techniques I would like to hear a little brief on that subject. Which are the possible techniques, computational power needed and so on.

If you have information on which part of the protein are interacting, and its sequence, you can search for homologous or similar sequences in other proteins. This would give you a list of potential targets and probably a similarity score which you can also use. Then you would have something to run simulations on.
And I think that with 3D structure you should be able to run simulations in silico in this case (eventhough it is probably more common to look for the ligands to the protein, not proteins to the ligand).
In silico simulations are computationally intensive, but I´ve also seen some being run on a laptop, so you don´t necessarily need access to HPC-clusters, or supercomputer.

This sounds like something you could possibly do, but is it safe to assume that the evolution of interactions is way faster than the sequence evolution? Is there any described method or paper using the approach you suggested?

@Kamil S Jaron And why we want to assume that?

You can find information in uniprot.org database (on structure, interactions, domains and motifs and similarity) and it also have supporting data on the data sources and literature citations.


text:
Integration of different microarray dataset to run GSEA
I'm planning to run a GSE Analysis on some microarray datasets. Until now what I have done is to pre-processed them using frma algorithm.
The reason I chose frma instead of gcrma or rma is that I wanted to merge these different dataset into one.
Here is how I'm planing to proceed:
i.e
dataset1 : 3x Controls , 3x Treatment1
dataset2 : 2x Controls , 2x Treatment2
dataset3 : 2x Controls , 2x Treatment3

After merging
7x Controls , 3x Treatment1 , 2x Treatment2 , 2x Treatment3

Then I will execute the GSEA for controls , Treatmen1 , Treatment2 and Treatment3 separably.
So the question is if it's more reasonable to follow the path that I described or to run the GSEA on the separate datasets (for example for each control ) and then get the intersect of the gene sets of these three repetitions.

Is there any reason you don't want to or can't do both?

What have you done to control that each dataset is (presumably) done under different conditions and thus have batch effects?  Are the controls the same in each dataset? Could you point me what is the frma algorithm ? It is the first time I read about it.

@Mitra : No there is no such a reason. I'm just asking to see if something is more reasonable or makes more statistical sense than the other.

@Llopis : At the merging stage I will remove the batch effects with comBat method (I could also try others as well). Here is a documentation for frma

Thanks, take into account that this will remove too the effect of each treatment. as it cannot distinguish between the batch effect and the treatment effect.

It makes more sense to evaluate by separate the pathway or gene set you want and see if in the three datasets result in a coherent message than to merge these datasets, as you will mix batch effects and treatment effects.
When in several datasets one looks for the same statistic it is usually considered meta analysis, not integration.


text:
Link Prediction in Bipartite Networks between biological pathways and drugs
Currently I'm building a big matrix (using microarray, mass spectrometry, RNAseq data) that consist by pathways (rows) and treatment/drugs (columns).
The values of that matrix are scores that describe the effect/relation of that specific treatment/drug on each pathway.
Now, some of the drugs is known to target specific disease (target disease) and thus is expected to have a high score with that disease's pathways. (A)
The rest of the drugs is known to target other diseases.(B)
My thought now is that, if I label the known drugs as "1" (or whatever), could I use a prediction model to find out new relationships (links) between drugs from the B class and the significant disease pathway from A ?
The whole thing, we can represent it as a bipartite network shown below.

For example, here the denoted as green, drugs are the known drugs belong to the class A that mentioned before and they are connected highly with some pathways of the target disease.
From the other hand blue and red ones are drugs taken from other cases (class B) that are used to target other diseases, but as my plot shows, blue seems to has an "affinity" to some of the pathways of the target disease and thus should be an alternative candidate for that disease. The red from the other side, doesn't have any "affinity" to any significant pathway and thus we can not consider it as possible candidate.
To find a solution to that kind of problem, I looked for algorithms and methods developed for social networks (e.g recommendation etc.) and mostly I read about random walks with return.
But before to proceed, I would like to know if someone else tried such an approach here and having more experience.

I am not sure if I understood correctly how you classify  the drugs. But what you attempt to do is similar to what they do here but in that article they use the drugs targeted to a molecule/pathway to find combination of drugs that are better for the disease (not an alternative drug), see the first image:

In that article they mention the pathway cross-talk inhibition (PCI), you might want to calculate it for your pathways to assess the relative reduction of the pathway efficiency.


text:
Interpreting Intergrative Genomic Viewer (IGV)
I was following a tutorial on "Tuxedo Genome Guided Transcriptome Assembly Workshop" and was wondering how to interpret the following:

From what I understand from 'Color Legends', the color blue represents something that is below normal and red represents something that is above normal.
Also, my understanding is that merged.gtf is the reconstructed transcript and genes.bed is the reference. They look to be similar for the most part, but seems to have breaks between. What's the appropriate way of interpreting this result?

You've apparently colored your alignments by read strand. In this case, red indicates "+ (watson) strand" and blue indicates "- (strand) strand". This strand association is determined by the orientation of read #1 in a pair (or just "the read" if you have single-end data). This isn't actually documented in the user manual, but you can find one of the Broad institute software engineers replying here.
Your next question is likely why the color of the reads and the orientation of the genes are opposite of each other. This is because in modern directional protocols, read #2 in a pair matches the orientation of the gene from which it arose.

Thank you for your reply. I'm not sure I understand what you're referring to by read #1 and read #2 based on the picture I attached.

Ignore that part if you don't understand it. The direct answer to your question is that the colors represent the strand from which the reads arose.


text:
Differences between NetAffx Hg-U133 Plus 2.0 Annotation file versions
From what I was able to observe, the annotation files (CSV) for NetAffx Hg-U133 Plus 2.0 is updated regularly.
However, from glancing at the CSV files (NA35: 2015 & NA22: 2007), I can't seem to notice any glaring difference between the two.
Does anyone know what information is actually being updated? I assume it may be gene identifiers for different databases, but wanted to make sure.
Thanks.
Below are the files I downloaded and examined from Affymetrix website (you need an account to download archived annotation files)

HG-U133_Plus_2 Annotations, CSV format, Release 22 (20 MB, 5/31/07)
HG-U133_Plus_2 Annotations, CSV format, Release 35 (33 MB, 4/16/15)

After checking the files myself, it does seem like gene identifiers indeed have been updated. Some that were missing have been added with identifiers and some have additional identifiers now.
Among some other things that were updated from 2007 include:

Genome version NCBI Version 36.1 => GRCh37
Alignments
Gene title
Gene symbol

Thought this would be a good reference for those who may be looking into older researches involving microarrays. I will be sure to check the files more carefully before posting questions like this one.

Nice you could answer yourself the question! (you can mark it as solved yourself) But I would still use the most up to date annotation of microarrays when using "old" data. However, I agree that some differences between analysis might be due to using different annotations.


text:
Highlight CpG sites with Texshade
I am using Texshade to produce figures for MSA. Upon reading the documentation I could not find anyway to highlight how conserved CpG sites within the alignment. Does anyone know if this is possible? If not is there another MSA visualizer I can use that can produce a figure like this. Thanks!

The best solution I could find was to use the \shaderegion function on all the sequences so it ended up looking something like this:
\shaderegion{1}{CG}{Red}{OrangeRed}

\shaderegion{2}{CG}{Red}{OrangeRed}

...

Then I cast the rest of shading similarity to grey so the CG pairs were more easily distinguishable.
\shadingcolors{grays}

Two examples one with \fingerprints and one regular.


text:
Major variability in speed of BLAST between proteins
I'm running blast+ searches for a list of proteins using specific parameters.
This is my search
 blastp -query protein.fasta -remote -db refseq_protein -out
protein_refseq_blast.txt -outfmt '6 qseqid sseqid qstart qend length evalue' -
entrez_query 'txid2 [ORGN] or txid4751 [ORGN]' -evalue 200000 -max_target_seqs 10000

For some proteins, this runs moderately quickly, ~5-10 mins. For other proteins, this command never stops running, even after an hour or so. I'm wondering if this variability in speed is due to any part of this choice in parameters, is simply variability in speed for the BLAST servers. Is it common for BLAST searches to take this long? Is there a way to run all of them in parallel using these parameters if this is the case?
Thanks.

FWIW, I remember experiencing very long waiting times when running remote blast from Biopython.

You have your evalue set incredibly high -evalue 200000, which might affect the time it takes to return a query.

The variability could be simply because the proteins are very different sizes but even 5-10 mins is very long for a single protein, as fugu says, your evalue is crazy! Probably reducing the evalue to a 'normal value' (<1) will resolve this

Try changing you evalue threshold to something like 0.001.
If you blasted a sequence against a database containing 10,000 sequences and got a hit with an evalue of 1, this would mean that, given the size of your database and length of protein, you would expect to find a match as good as (or better than) the hit returned once purely by chance. If the evalue threshold in your query was set to 0.01 (for example) then this hit would not be returned.
If you got a hit with an evalue of 200,000 this would mean you would expect to see 200,000 matches with a similar score simply by chance. If your evalue threshold were set to 200000, then all of the hits with an evalue <=200000 would be considered significant. If you set it to 0.001 then only those hits with evalues <= 0.001 would be considered significant.
This means that (depending on the size of the data base) almost every sequence in your query set will be considered as significant given your threshold, which will slow down your search hugely.
See the the section on E values here

The evalue associated with a blast HSP (High-scoring Segment Pair) is a measure of how often you would expect to find such an HSP in a database of this size by pure chance. For instance, if you blast the sequence ACTG against the entirety of the NCBI's nr database, you would expect several hundred thousand hits by pure chance. That is, several hundred thousand sequences that can be aligned to your query but are in no way related to it (not actual homologs). Those hits don't actually mean anything, so you use the evalue to limit the results to those that are unlikely to be found by chance and, therefore, more likely to have biological meaning if hits are found.
So, setting the evalue to 200000 means you are interested even in HSPs whose characteristics (score, length, etc) mean you would expect to find <=200000 HSPs with such characteristics by pure chance. This evalue is pointless.
Because your e-value is so high, BLAST will report absolutely anything and everything it finds, and the vast majority of your results will be meaningless and irrelevant. If you want a permissive e-value, you could use 1 or even 10 but 200000 is way too much. You usually only want negative e-values, cases where you wouldn't expect to find even one such HSP by chance. For highly dissimilar sequences, I would recommend you use 0.1 or something of that scale.
Changing the evalue will speed things up considerably. That said, you should always expect different times for different queries since the time it takes to find hits will always depend on the number of hits that can be found. If your sequence can be aligned to many of the sequences in the database, it will take longer to fetch all results.


text:
how to set database other than nr for remote blast+ search
I am attempting to run a BLAST search remotely using BLAST+.
I can get search to work correctly at the command line with the following commands:
blastp -query proteins.fasta -remote -db nr -out proteins_nr.txt -outfmt 6 -evalue 1e-30

However, I would like to us the remote database titled "Microbial proteins from nr" which is what would be used for the microbial blast search accessing BLAST from the NCBI site.
I'm struggling to find the correct code to access this database instead of just "nr".
Is there a list with the codes of different databases?
Thanks in advance

I believe you're looking for env_nr? It's listed as such, under Metagenomic proteins in the blastp webpage. It appears that the word within brackets should be supplied alongside the -db parameter. A quick test with a dummy amino acid fasta file does turn up a result to a valid NCBI protein accession.

Edit: I followed up on this with a little more digging around, and it appears the database you're looking for is called Microbial_proteins. I ended up quote-searching Microbial proteins from nr on Google and exactly 2 hits turn up. The first hit, I think, is what the OP posted a screenshot of in a comment to this answer. The second is a list of databases along with keywords from a Spain-based bioinformatics firm: http://data.biobam.com/ncbi_blast_dbs_protein.pdf. I retried the dummy example with the Microbial_proteins database, and it seems to work. Again, I am not sure why I can't easily find this on an NCBI resource. Also, I am curious as to how different this database is from env_nr.

@flatley176 thank you for response. That answers what I can specify with -db. Unfortunately I'm not looking to blast against env_nr. I would like to blast against the website database using blast+ and specify this database:

is there a way to specify this for a -remote blast using blast+
thank you.

I do not know how Microbial proteins from nr differs from the env_nr database. Page 3 of the BLAST manual (ftp://ftp.ncbi.nlm.nih.gov/pub/factsheets/HowTo_BLASTGuide.pdf) lists the available databases - the only reference to the word "microbial" is to the 16s rRNA. That said, I dug around a little more, and it appears there IS a solution - the database is called Microbial_proteins. Will update my answer.


text:
Does NCBI's blast API block my IP?
I am trying to run a blatn command from both my laptop and within a google colab notebook. I am not sure why, but this command runs properly once and then, on the second execution, it gets stuck for hours without a reply.
Is it possible that NCBI blocks my IP? I've tried also to run web_blast.pl and this didn't work either.
biopython:
from Bio.Blast import NCBIWWW

fasta = '>seq0\n'+'AACCACCTGACGGGCTCTGT'
blast_res = NCBIWWW.qblast('blastn', 'nt', fasta)
print(blast_res.read())

web_blast.pl command:
web_blast.pl blastn nt input.fa > blast.out
Note, if this is due to throttling by NCBI, I am not sure how do they expect people to query their web service if they block you after one call.

If you have a high usage department might explain it. Try via your own ip

You should check the User Guide:

The NCBI BLAST servers are a shared resource. We give priority to
interactive users. In order to ensure availability of the service to
the entire community, we may limit searches for some high volume
users. Interactive users of the NCBI webpages through a web browser
should not encounter problems. We will move searches of users who
submit more than 100 searches in a 24 hour period to a slower queue,
or, in extreme cases, will block the requests. To avoid problems, API
users should comply with the following guidelines:

Do not contact the server more often than once every 10 seconds.

Do not poll for any single RID more often than once a minute.

Use the URL parameter email and tool, so that the NCBI can contact you
if there is a problem.

Run scripts weekends or between 9 pm and 5 am Eastern time on weekdays
if more than 50 searches will be submitted.


text:
How to get Uniprot and swissprot cross references using Biopython
I am struggling to use biopython to gather cross references such as GO annotation from Uniprot/swissprot text files ("DR" in the .txt files) using biopython.
From the source code:
def _read_dr(record, value):
    cols = value.rstrip(".").split('; ')
    record.cross_references.append(tuple(cols))

This leads me to think that the following code should work, as this is how other features like annotation and features work.
for record in SeqIO.parse(filename, "swiss"):
    print(record.cross_references) #does not work, should return a list of DR rows
    print(record.features) #works, returns a list of parsed FT rows

I believe the attribute you should be looking for is dbxrefs:
[x for x in record.dbxrefs if x.startswith('GO:')]

Apparently SeqIO handles uniprot text files differently than the SwissProt module:
from Bio import SwissProt
record = SwissProt.read(open('single_record_from_uniprot.txt'))
record.cross_references

So you should change you code to use the SwissProt module to get the behavior you expect


text:
Differential Expression With 2 Treatments
I have a expression data from a small cohort of samples taken at baseline and after 2 independent treatments. I can do differential expression contrasting T1 and T2 or I can contrast T1 vs baseline and T2 vs baseline and look at the differences. What is the difference between these two analyses? Are they both valid? If so what inferences can be drawn from each?

If you're interested in looking at the differences between two treatments then you'll end up wanting to do both a direct contrast as well as the individual comparisons to baseline.
The direct contrast will give you the genes actually differentially expressed between the two conditions. In practice, you may want to filter this a bit so you only have genes differentially expressed vs. baseline in at least one condition (e.g., to get rid of genes only slightly higher due to T1 and slightly lower in T2, but not different enough in either case to be DE). Use a pretty lax p-value threshold (e.g., 0.1 or 0.2) for this.

Thanks. What is your opinion on genes that may be significant in just 1 of the comparisons to baseline but not significant in the direct T1 vs T2 contrast?

It there's no difference in the T1 vs. T2 contrast then there's no difference. It doesn't matter if you happen to get a significant difference in one of the conditions vs. baseline. You only want to use the comparisons to baseline for filtering your T1 vs. T2 results (assuming you have enough DE genes to bother doing so).

One more question. So if I wanted to do a pathway analysis and compare the two treatments that way. When I generate my gene lists how should I treat genes that are DE in only one treatment vs baseline but not DE in the direct T1 vs T2 contrast?

Treat them as not DE.

So any genes that have similar responses from baseline in both treatments get discarded? Wouldn't you be missing a lot of context if a gene has to be DE vs baseline in both treatments and also DE directly with each other to make either final list?

Yes, any genes with similar responses should get discarded, since they're not different. A gene doesn't have to be DE in either treatment vs. baseline to be DE between them. That's also why I suggested using a very lax p-value threshold if you want to do any filtering of the T2 vs. T1 results.

Sorry, I meant doing a pathway analysis of each treatment vs baseline separately to find the raw response of each treatment. Since they are different from baseline I would still include them in only their respective pathway analysis?

Yes, certainly include them in such cases.

Cool, thanks a lot for your help. I really appreciate it.

Depending upon the type of treatment used the set of DEGs will change.If the treatments have similar kind of effect you will get a small list (less variable genes will have higher p-val) using a cutoff of p<0.05. So, it's better to start with control vs treated comparisons then T1 vs T2.
Compare the lists form CTRL vs T1 and CTRL vs T2 you will get the genes that are expressed in both conditions as well as unique to individual treatments.

If a gene has a p-value of 0.049 in CTRL vs. T1 and 0.05 in CTRL vs. T2 you would be categorizing it as uniquely DE in one condition, when there's no actual difference between the conditions.

Yes, it's a tricky situation.

I agree it is a trick situation in my case it is modes of exercise so they are extremely similar responses. The T1 vs T2 comparison yields almost nothing where the vs baseline comparisons are quite different.


text:
Appropriate tool or algorithm for sloppy alignment of degenerate bases
I have an optimization problem where I have a degenerate nucleotide sequence I want to align against subsets of a reference genome (exons, specifically, to make the problem more tractable).
The sequence itself contains subregions of variable length.
The sequence looks like this, containing degenerate bases:
           L           C           R
5' - SSS......SSS|NNN.....NNN|WWW......WWW - 3'

Per IUPAC notation, S denotes C or G, N denotes any base, and W denotes A or T.
I am marking each block L, C, and R to broadly denote "left", "center", and "right" regions, oriented from 5' to 3' ("forward"-strand oriented).
The widths of the L and R blocks can be from 20 to 25 bases long. The width of the C block can be 14 to 16 bases long.
One programmatic way I could think to do this would be to brute-force a solution, generating all combinations of sequences of varied lengths, blastn-aligning each of such combinations to the sequence for the exon, and then looking for the best alignments among such sequences.
The output would be the scores for regions from sequences that align, to be able to say, from point X to point Y of the genome, a sequence here aligns better than other candidate sequences.
To avoid reinventing tools, and avoid brute-force methods, I'd like to know if there is a tool or even a wrapper around tools, which does this, in case it is already written, or another approach that is more elegant/efficient.

I believe mafft can align degenerate sequences.  From the documentation for version 7 at: https://mafft.cbrc.jp/alignment/software/anysymbol.html
Acceptable nucleotide symbols in the default mode:
a, c, g, t, u
Ambiguous nucleotides (r, y, w, s, k, m, d, v, h, b; IUPAC-IUB codes) can be used and are scored as:
score(r,a) = ( score(a,a) + score(g,a) ) / 2
score(r,g) = ( score(a,g) + score(g,g) ) / 2
score(r,t) = ( score(a,t) + score(g,t) ) / 2
score(r,c) = ( score(a,c) + score(g,c) ) / 2
...
score(r,y) = ( score(r,t) + score(r,c) ) / 2
...
score(r,r) = ( score(a,a) + score(g,g) ) / 2
score(y,y) = ( score(c,c) + score(t,t) ) / 2

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

Can you expand your answer to include some support material?

I thought that a nucleotide blast could do the same. Just wondering if these tools also support querying against variable-width subsequences, or if wrappers around such tools can do this more efficiently in generalized cases. I have worked on parallelized analysis of combinations of substrings, in the meantime. It is brute-force but the scope of the problem seems tractable.

I'm not sure I quite understand your point about variable-width subsequences - as far as I understand, these tools tend to find the longest alignable subsequences they can from the two input sequences, so I guess if you want to test lost of different widths a wrapper would be the way to go.

@AlexReynolds I don't see why this isn't a 'normal' alignment problem where you could use blast or similar but with relaxed parameters (maybe permissive gap penalties to handle "variable-width" subsequences etc.). Perhaps a little more biological context would help us understand better too


text:
Efficiently concatenating disjoint, ordered bigBed files
Given a set of bigBed files, each of which are disjoint to each other and in lexicographical ("sort-bed") sort order, I would like to know if there is a way to concatenate these archives efficiently.
Presuming I know the order in which to specify files (to constrain the order of intervals), is there a way to concatenate them directly, i.e., without first extracting them to bed, concatenating the extracted files, and remaking a bigBed file?

No, there's no simple way to do that without someone writing a program tailored toward it that handles the footer section and recreates the zoom levels that follow the compressed blocks of intervals. If you wanted to write such a program, modifying bigWigCat from Kent's tools would be a reasonable tack.

Seems reasonable, thanks.


text:
Error in checkFullRank(modelMatrix) deseq2
Trying to use deseq2 for differential expression analysis (rna-seq) between three groups and also account for batch effect as the control were sequenced at a different time point.

control: con
sample with mutation A: mutA
sample with mutation B: mutB

here is my design
   sample condition batch
1     C_1         con     1
2     C_2         con     1
3     C_5         con     1
4     C_3         con     1
5     C_4         con     1
6     M_6        mutA     2
7     M_2        mutA     2
8     M_5        mutA     2
9     M_1        mutA     2
10    M_4        mutA     2
11    M_3        mutA     2
12   MA_2        mutB     2
13   MA_6        mutB     2
14   MA_5        mutB     2
15   MA_4        mutB     2
16   MA_3        mutB     2
17   MA_1        mutB     2

when i create matrix:
dds <- DESeqDataSetFromMatrix(countData = data, colData = coldata, design = ~batch + condition)

I see this error:
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

I also tried to compare only two conditions subsetting the original data, i still see the same error what am i missing?

can you please post the format which ran successfully in DESeq. I have a similar issue

See the second example from the vignette, which is exactly your situation. The software cannot differentiate between the batch effects vs. condition effects because the groups are exactly the same.

cannnot re do the experiment, use the nested column is it?

elegantly explained just had to a little more search, can i delete this question?

Sorry, I don't understand your question in the first comment. Sure, you can delete it if you'd like.

no the link you posted has accurate answer, i had to do a little more research before posting a question how should i design was my question


text:
ribo interval list for picard ensembl mm9 gtf
Could some one please help me with understanding how to generate the ribosomal interval list that is required to when using picard metrics this step
java -jar picard.jar CollectRnaSeqMetrics \
      I=input.bam \
      O=output.RNA_Metrics \
      REF_FLAT=ref_flat.txt \
      STRAND=SECOND_READ_TRANSCRIPTION_STRAND \
      RIBOSOMAL_INTERVALS=ribosomal.interval_list

I have a couple of bam files that were aligned using mm9 ensembl gtf as annotation, i did a lot of search and found this link Ribosomal Intervals For Collectrnaseqmetrics but this link or any other source does not explain how to produce or find ribo interval list from ensembl annotation

At least with Ensembl annotations, you can grep rRNA foo.gtf > rRNA.gtf to get the annotated rRNAs. Note however, that the 45S cassette isn't in the mouse reference genome, so you'll never have alignments to 18S, 5.8S, or 28S rRNA. I always get the feeling that the Collect*Metrics programs are made with human sequences in mind, since the cassette is included in the human reference genome.
Anyway, while you can use the annotation and reference sequence to make picard happy, make sure to append the Rn45S sequence to your genome if you want to actually measure rRNA levels. In my experience this is only useful for rRNA-depleted samples, since the rates can vary wildly there.
Edit: To make a file compatible with picard's "interval list", one can do the following with the GTF file:
grep rRNA foo.gtf | awk 'BEGIN{OFS="\t"}{print $1, $4-1, $5}' > foo.bed

According to the picard documentation a BED file can be used as input and the above awk code extracts the appropriate columns and switches to 0-based coordinates used in the BED format.

Interesting, can I directly use the output of grep rRNA foo.gtf > rRNA.gtf as input file? or is there any further modification necessary? I did not understand when you say make picard happy I am using mm9 reference and yes the sequencing done was using ribo depleted technique (could you please point towards some literature or steps)

@novicebioinforesearcher I've updated my reply with the conversion to BED, which allegedly should work. "make picard happy" is another way of saying "give picard a file it can use as input". The results would still not be terribly useful without adding in the Rn45S sequence.

I never came across Rn45S concept (thanks for pointing it out) and have not seen any literature in the tools that use such list what would be your thoughts as to why they would not be speak about it?

@novicebioinforesearcher Many reference genomes (in particular, human and fly) contain it already, so there'd be no need to use Rn45S there. I don't know why it was excluded from the mouse genome, it's one of the things that's always bothered me.


text:
Difference between CPM and TPM and which one for downstream analysis?
What the difference between TPM and CPM when dealing with RNA seq data?
What metrics would you use if you have to perform some down stream analysis other than Differential expression for eg.
Clustering analysis using Hclust function and then plotting heat map to find differences in terms of expression levels, correlation and pca
Is it wrong to use TPM for such analysis, if yes then when does one use TPM versus CPM.

You can find the various equations in this oft-cited blog post from Harold Pimentel. CPM is basically depth-normalized counts, whereas TPM is length-normalized (and then normalized by the length-normalized values of the other genes).
If one has to choose between those two choices one typically chooses TPM for most things, since generally the length normalization is handy. Realistically, you probably want log(TPM) since otherwise noise in your most highly expressed genes dominates over small expression signals.

in a given case if one would trim adapters from paired end sequenced rna seq data, this would result in different read length, as you say above TPM is length normalized does it mean that this difference in read length is taken into consideration?

@novicebioinforesearcher No, transcript/gene length, which correlates (somewhat) with counts and therefore will tend to drive clustering unless handled in a reasonable way.

Rob Patro also wrote a pretty good article about the topic as well: http://robpatro.com/blog/?p=235

[link] (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8220791/)

Neither CPM nor TPM are well suited here, because neither performs robust cross-sample normalisation (see the blog post Devon linked to).
DESeq2 provides two robust log-space normalisation methods for downstream analysis, the regularised log (rlog), and the variance stabilising transformation (vst). The DESeq2 vignette explains how to use these for things like hclust.

On a more general note, CPM does not account for transcript length differences, while TPM does. If the choice is between TPM and CPM I would therefore use TPM. However, if you are only comparing the same transcripts across experiments, the transcript length is actually invariant so it doesn’t matter (but CPM is still not a good cross-experiment normalisation).

I am confused now , so what is the use of TPM, why does one produce it and when or where do you use it? , in other words which tools, analysis in RNA seq would you use TPM if everything revolves around using counts and pushing it through DESeq2

Tools produce TPMs because they don’t have the information (= the other samples) necessary to perform cross-sample normalisation. Lacking that, TPM is the best they can do. TPM is also useful for within-sample comparisons: It can give you an accurate estimate of how much genes are expressed in a given sample relative to each other.

Plz Correct me if am wrong here, given an experimental design different cell types from a normal mouse say 4 cell types (3 replicates each), sequenced using same library prep but may be at different times. Aim would be to check for set of cell type specific transcripts you would use TPM where as if you want to add a significance parameter (need a pvalue) one would use raw counts based analysis?  I guess the confusion for me is when we use the word  "expression" . People use TPM and call it expression, also use raw counts and call it expression

Both are estimates for expression, given the data. Your use-case sounds reasonable although I would generally prefer to determine “cell type specific transcripts” by comparing different cell types, rather than based solely on a single sample. Which would imply performing differential expression analysis.

Which would imply performing differential expression analysis on?

@novicebioinforesearcher On whatever data sets you wish to compare. It rarely (if ever!) makes sense to describe a gene as being cell type specific without saying “as opposed to these other cell types”. For example, a cell type specific gene may nevertheless be lowly expressed: as long as it’s entirely absent in other cells, it’s cell type specific. This is in fact often the case. You therefore cannot characterise many cell type specific genes without comparing different cell types.

Is VST/RLOG normalised for transcript length? If not is it possible to get length normalised VST/RLOG?

@rmf No, they do not normalise for transcript length; both functions merely change the distribution shape of the counts into something closer to linear. For those purposes where you’d use rlog/vst, accounting for transcript length, is normally not important. You could however apply further transformation (let’s call it rlog-TMP), if you have an application where you need cross-sample as well as within-sample normalised values.

Would you say it's ok to do something like (vst/length)*(10^6) and use that for heatmaps where I actually want to compare expression of one gene to another gene.

@rmf Yes but in a heatmap you usually scale by (gene) row anyway (either explicitly or the plotting function performs the division internally) so the normalisation per transcript length will be strictly a no-op.


text:
How to count reads in bam per bed interval with bedtools
I recently installed Ubuntu 16.04 (because I was still using 12.04). But it seems my bedtools scripts don't work properly anymore. I can't figure out how to use the new bedtools for my old ways. What I want to do is get the number of reads from a bam file, per interval from a bed file. It was very simple with my old version of Ubuntu and bedtools:
bedtools coverage -abam file.bam -b All_peaks.bed > file.cov.txt

But this doesn't seem to work anymore. I used to get small files as results (~3MB), but now huge files are created (~4GB) including read names??
I am using bedtools v2.25 now.

The order of -a and -b switched at some point. You want:
bedtools coverage -a All_peaks.bed -b file.bam > file.cov.txt

For reference, this is the end of the help output in version 2.25:
Default Output:
     After each entry in A, reports:
       1) The number of features in B that overlapped the A interval.
       2) The number of bases in A that had non-zero coverage.
       3) The length of the entry in A.
       4) The fraction of bases in A that had non-zero coverage.

And this is the equivalent output from version 2.19:
Default Output:
     After each entry in B, reports:
       1) The number of features in A that overlapped the B interval.
       2) The number of bases in B that had non-zero coverage.
       3) The length of the entry in B.
       4) The fraction of bases in B that had non-zero coverage.

Update: The change in behavior happened in version 2.24:

We have changed the behavior of the coverage tool such that it is consistent with the other tools. Specifically, coverage is now
  computed for the intervals in the A file based on the overlaps with the B file, rather than vice versa.

Thanks Devon! I wasn't aware of the 'switch'. Switching -a and -b solves indeed my problem!


text:
install bowtie2 from sources cannot find -ltbb
I am trying to install bowtie2 aligner from sources without the root access.
Bowtie2 needs tbb package to be installed and it is recommended to install oneTBB.
# install one TBB
mkdir -p SHOME/soft
cd $HOME/soft
git clone https://github.com/oneapi-src/oneTBB.git
cd oneTBB
make

Then I add C_INCLUDE_PATH to .bashrc:
echo "export C_INCLUDE_PATH=$HOME/soft/oneTBB/include" >> ~./bashrc

and continue to install bowtie2:
# install bowtie2
mkdir -p $HOME/soft/bowtie/source
cd $HOME/soft/bowtie/source
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.1/bowtie2-2.4.1-source.zip
unzip bowtie2-2.4.1-source.zip -d ..
cd ../bowtie2-2.4.1
make

TBB library is recognized but in the end I get this error:
/usr/bin/ld: cannot find -ltbb
/usr/bin/ld: cannot find -ltbbmalloc_proxy
collect2: error: ld returned 1 exit status

ltbb and ltbballoc_proxy are not in TBB but libtbb and libtbballoc_proxy instead.
I also tryed to add variables to .bashrc:
TBB_LIB_PATH=$HOME/soft/oneTBB
TBB_INCLUDE_PATH=$HOME/soft/oneTBB/include
TBB_RELEASE_DIR=$HOME/soft/oneTBB/build/home/suvar/apps/oneTBB/build/linux_intel64_gcc_cc9.3.0_libc2.31_kernel5.4.0_release

but seems it does not help.
I need help to complete the installation of bowtie2 on our server.

Did you remember to source your bashrc before compiling?

sure: source ~./bashrc

Do you really need to compile Bowtie2? There are compiled binaries for Linux and MaxOsX available:
https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.1/bowtie2-2.4.1-linux-x86_64.zip
https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.1/bowtie2-2.4.1-macos-x86_64.zip

If you are using Ubuntu you can install tbb with sudo apt install libtbb-dev, in MacOS you can install it with Homebrew brew install tbb.

You need to give the path to the source files for the compiler and the full path to the compiled libraries (all parent directories included) for the linker.  You also need to use either CPATH or CPLUS_INCLUDE_PATH instead of C_INCLUDE_PATH to point to the oneTBB source, since bowtie2 is C++ and C_INCLUDE_PATH will only have an effect for C (see Environment Variables Affecting GCC in the GCC documentation for more info.)
For example like this, inside bowtie2-2.4.1, without making any lasting changes to the shell's setup in this case:
$ ONETBB=$HOME/soft/oneTBB
$ export LIBRARY_PATH=$ONETBB/build/linux_intel64_gcc_cc5.4.0_libc2.23_kernel4.4.0_release
$ export LD_LIBRARY_PATH=$LIBRARY_PATH
$ export CPATH=$ONETBB/include
$ make

After it's built you'll still need the library path to oneTBB defined (via LD_LIBRARY_PATH) so bowtie2 can find it at run-time.
But if you don't have a specific reason to compile from source I'd recommend avoiding it.  You could go with what h.mon suggested and use an existing binary, or use a package manger to install it like zorbax suggested.  Conda has a package for bowtie2 in the bioconda channel, up through and beyond version 2.4.1, and you can use that as a package manager without root access.


text:
ncbi...eutils/esearch.fcgi: How to search now in the snp database?
According to the 2017(!) textbook, it worked at that time, for instance,
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=autism
At this time, this does not work and for each conceivable term it will find 0 records.
What is the right way to search in snp, valid in December 2018?

Welcome to the site Slepecky. Could you provide a reference to the textbook you mention? Maybe there is an updated version someone can look up. Have you tried looking the eutils manuals for this task?

Also, please explain what you are trying to do so we can come up with alternatives.

eUtilities for dbSNP is in the process of being retired. See this post from NCBI. You may want to check out the new API.
From the NCBI Insights blog post:

Entrez search is available only for human. We will retire all Eutils report formats (RS docsum XML and ASN.1, FASTA, and FLAT text) later this year. Only eSummary XML report will be supported.
Batch Query Service for all organisms will be retired October 2018. Users of human dbSNP data should transition to the new API.
We will place non-human FTP data in the archive directory. Please see the phase-out plan announced last year.


text:
Is there a database of protein sequences/structure along with their melting temperature?
Such databases have been constructed for example in here or here. But I can't find them anywhere online. Can someone point me to where I can find one available online?
Otherwise I'll have to write to these people and I don't know when/whether I'll get a reply.

ProTherm2 contains thermodynamic properties of proteins and mutants. https://www.iitm.ac.in/bioinfo/ProTherm2/ .
NAR Database issue advertises it as containing melting temperatures.
For the reference - found it with https://biokeanos.com/search?q=protein+%2Bmelting

Thanks! What amuses me is that back when I needed the database, we did contacted the authors of ProTherm for the data and used it. Seems like it's the only available database for this purpose, right now.


text:
How do I find a list of PDB structures bound to a peptide inhibitor or peptide-like inhibitor?
I looked up existing databases but ones like iPPI-DB or enter link description here either provide non-peptide inhibitors or only sequences. I need structures of proteins bound to other proteins, the latter of which is an inhibitor of the former, and the complex has a structure deposited in the PDB database. Here's an example: PDB:3U1J
In case there is no existing database in literature, it'd be great even if I can screen for such PDB entries using some code, but I can't come up with an algorithm for it. Thank you!

Welcome to the site @simple_dozer! I'm not sure if there is a code name to know if a protein is an inhibitor, but you can start by subsetting all the structures with at least two genes

The quickest way to screen the PDB database is Advanced Search from either RCSB or PDBe.
Since you need entries with 2+ proteins you can set the number of protein chains to at least 2 in the RCSB search.  Or you can select "protein/protein complex" in Assembly composition in the PDBe search.
It's less clear how to search for inhibitors. If you'd look for the "inhibitor" word in the molecule names you'd get (together with the previous criterium) 1200+ entries.


text:
What methods should I use from PythonCyc API to query metabolites in BioCyc database?
I am using PythonCyc API in order to write a query for metabolites in BioCyc. The purpose of this API is to communicate with the database software of BioCyc- Pathway Tools. Pathway Tools is in lisp therefore, PythonCyc creates a bridge between python and common lisp. To use this this API one must first create a PGDB object with a specified organism ID (orgid). In the example below I create a PGDB with orgid "meta". After this, I am able to call methods from PythonCyc with this object:
import sys
import pythoncyc
#this creates PGDB object associated with meta(MetaCyc)
meta = pythoncyc.select_organism('meta')
#lists pathways of specified compound
pathways = meta.pathways_of_compound('sucrose')
print pathways

The query for the metabolite above, 'sucrose', provides a list a pathways:
[u'|PWY-7347|', u'|SUCSYN-PWY|', u'|PWY-7238|', u'|PWY-6527|', u'|PWY-5343|', u'|PWY-3801|', u'|PWY-7345|', u'|PWY-862|', u'|SUCUTIL-PWY|', u'|PWY66-373|', u'|PWY-621|', u'|PWY-822|', u'|SUCROSEUTIL2-PWY|', u'|PWY-6525|', u'|PWY-6524|', u'|PWY-5337|', u'|PWY-5384|']

However if I switch this metabolite name to a common amino acid name, such as 'valine':
import sys
import pythoncyc
#this creates PGDB object associated with meta(MetaCyc)
meta = pythoncyc.select_organism('meta')
#lists pathways of specified compound
pathways = meta.pathways_of_compound('valine')
print pathways

The query provides an error message that states "non coercible frame", which means it cannot find this ID. This is a common metabolite, just as sucrose, which should very well have pathway entries in database, however there are none found with this method. I have also tried synonyms of valine and other methods such as "all_pathways(cpds)" stated in API, which give me the same error message.
What is the reason for this? Does BioCyc not list any pathways for amino acids? Or am I using an incorrect method to access amino acid information?
Link to PythonCyc API

Hi astrika, thanks for your question and welcome to Bioinformatics Stack Exchange. I've had a glance over this question, and it looks like it has the potential to get good answers: it is well-researched, with a good story, and a specific question. I look forward to seeing more questions from you in the future.

Pathway-tools has a unique identifier (called a Frame-ID) for each metabolite in its database.  You can get the list of all compounds with this api call:
all_cpds = meta.get_class_all_instances('|Compounds|')
for each compound, you can also get its common name, or all its synonyms:

cpd_names = {}
for cpd_frame_id in all_cpds:
      cpd_names[cpd_frame_id] = meta.get_name_string(cpd_frame_id)

UPD: As suggested by Llopis, maybe the correct identifier is VAL?
I think there is no metabolite with such name in the database. I tried their search, the correct name seems to be L-valine. I expect that they use similar naming for other amino acids. If you cannot find something, you can use their compound search through the website.

I have tried 'L-valine' in the script above and it has also gave me the same error. I am aware it is possible to search in the site, however for my purpose I must write a query. My purpose is to obtain possible pathway networks for a given metabolics dataset

@astrika maybe the name is Protein-L-valine, or VAL. One thing is how they store things an another is how they show them. Sorry I didn't test how to search using python


text:
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.

Thanks Devon! multiBamSummary tool produces a scalingFactor table for respective BAMs. Does it make sense to first scale with these scaling factors, create bigwig using bamConverge and then use these bigwigs for comparison using bigwigcompare? I am wondering if scaled files go to bigwigcompare, will it make for a better comparison across comparisons.

Yes, that's another good way of approaching things. bigwigCompare does assume that the input bigWigs are scaled appropriately and the procedure you suggest will create a nicely robust comparison.

Thanks! You mentioned plotFingerprint for seeing good separation. I did the plotCorrelation heatmap. Is there some approximate threshold for "good separation"? Let's say I have a correlation between 0.5 and 0.7 across different samples. Will this make them ok for SES?

There's no substitute for plotFingerprint here.

Devon,
Can we use scalingFactor table (from multiBamSummary) to create new scaled BAMs instead of bigwig/bedgraph? I want to perform peak calling only on normalized/scaled BAMs so I can properly compare peaks/footprints across different samples. I know I can downsample using samtools but some scaling factors are greater than 1 so not sure how to "upscale" it.
Thanks again!

So  I found a post where you suggest not to upsample:
https://www.biostars.org/p/365825/
But lets say after running multiBamSummary, i get a table that looks like this:
A     1.1215
B     1.4657
C      1.3783
D      1.0529
E      0.6131
F    0.5919
G       0.9988
H   1.9739
I   1.0636

In order to make these samples comparable, do you have any suggestions on how do I go about downsampling everything without upsampling?

At least some peak callers can accept bedGraph files, which you can normalize. The issue with normalizing BAM files is that you're quite limited in terms of fidelity (e.g., you can't remove a tenth of a read from a region).


text:
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?

Which kind of physical interaction are you looking for: an hydrogen bound, a covalent bond, a ionic bond...?

There is no different for me which kind of physical interaction. The least distance and the farthest no different just physical

The CASP experiment uses an 8 Angstrom distance between C-beta atoms (C-alpha for glycine) to determine residue-residue contacts.

Use a distance cutoff of 12Å between Cᵦ atoms.

If this is part of your previous answer, consider editing it, otherwise it would be nice to provide bit more explanation.

Please could you clarify what types of interactions can happen at such distances, as I am currently unaware of interaction being possible at such a long distance?

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".

You should read my paper, where I did a comprehensive parameter sweep across different thresholds and definitions and benchmarked the resulting empirical potentials on a decoy ranking challenge:
The relevant phrase from the abstract is "compared using 90 different definitions of residue-residue contact". However, it boils down (as always) to the question, "what do you want to do".
Cite me if you use my paper for your research.

Most papers give 5 or 6 angstrom distance for physical interactions. I don't know this distance is between which atoms exactly

People use what other people use, which is why I did a 'parameter sweep'. I've seen other work using a strict 'information maximization' criteria, which was nice.


text:
Finding residues physically interacting in a PDB structure and sequence related to that structure
I am going to describle the problem completely. First I should say that there are two kinds of interactions: transient and stable interactions. We are supposed to work on stable interaction like interactions between two monomers in a heterodimer. In a heterodimer there are two chains in which there are some residues between monomers which take part in the interaction of monomers and building the structure of heterodimer and also there are other residues interacting inside each monomer for building the structure of that monomer. There are some PDB IDs related to heterodimers in PDB. Now imagine that we are trying to show which residues are physically interacting between two monomers in a heterodimer and then we are trying two show the site of these physically interacting residues in the sequence related to that special structure. For that, we need to find the sequence related to that structure and we can download it from PDB. But if we want to make a multiple sequence alignment of that sequence, first we need to find the refseq sequence of that and then we must blast the sequence and find homologues of that sequence. Here there are some challenges: first, we must find residues interacting in the structure of heterodimer between two monomers. Second, we must find the sequences of monomer chains of the heterodimer and map the interacting residues in the structure to the sequence. Third, we know that when a structure is solved and its monomers are sequenced, may be they are not able to completely sequence the heterodimer and residue numbers are maybe different from the related refseq sequence. I would like to know how can we find the related refseq sequence (or sequences, in the case of heterodimer) and then how can we map the physically interacting residues in the structure of PDB ID to refseq sequence (finding the sites of these interacting residues in that refseq sequence) regarding three challenges mentioned above. Imagine that I have 100 PDB IDs of 100 nonredundant heterodimers and I would like to find and download the structures, sequence and finally refseq sequence. What should I do?

Dear jaslibra. Thank you for your nice comment

My reputation isn't high enough to make a comment, but have you considered a modelling software such as PyMol or Chimera? You can execute a command line script with softwares such as these and these have the functionality to predict H-bonds.
In the PDB files there is information regarding the GenBank and UNIPROT info of the protein, and you can find more information here:
http://www.wwpdb.org/documentation/file-format-content/format33/sect3.html.
If you want RefSeq ID, you can use GenBank info you extract and convert by piping into a biomaRt (Bioconducor) script, although maybe there is a more direct way.
As for finding which residues in the sequence are interacting, I'm fairly certain that the software packages like PyMol or Chimera will give you information about the interacting residues.
From what I understand, you will want to filter out only the residues that interact and are on different chains (e.g, Chain A vs. Chain B). This process will require you to work with the syntax of the different software packages when you are scripting. For example, Chimera has a basic primer here:
https://www.cgl.ucsf.edu/chimera/current/docs/ProgrammersGuide/basicPrimer.html
Note that I only gave a reference to Chimera because I am more comfortable with that then other packages and there is likely a more efficient way to do this. With Chimera, you can find the position within the chain of residues participating in the interaction and then you'd store these in a file to do whatever you want later.
I hope this helps.

Hi jaslibra, thanks for your comment, and welcome to Bioinformatics Stack Exchange. This comment is well-structured, and looks long enough and verbose enough to be an answer for this question, so I haven't flagged it as a comment. I look forward to reading more comments and answers such as this in the future.


text:
Finding homologs of a protein sequence
I have a refseq ID of a protein from E.coli and I want to find homologs of this protein. I ran Blast against refseq database but I got a lot of sequences most of which were from Ecoli again. I decided to run PSI-Blast to get more divergent species, but I do not exactly know if my result are real homologs or false positives. what is your idea for finding homologous protein sequences from more divergent species? And what can I do for selecting the real (but not false positive) hits?

False positive can be judged by the E and P-values so remember to use them.

I would suggest that if you can robustly obtain the same topology using a variety of phylogenetic methods and models on your proteins, there is a chance that the homology is true. Unfortunately, this might be difficult to obtain when dealing with quite distant species, even if the homology is true.

It depends on what you're looking for. If you're just looking for sequence homology, then you can simply pick the best hits from a blast search. If, however, you are referring to functional homology, if you are looking for the protein which has the same functions as your query, then it's more complicated.
Sequence homology is not enough to infer functional homology. For example, you can have cases of gene duplication and subsequent functional divergence. Such paralogs are still homologs (paralogs are a subset of homologs), but they don't necessarily have the same function. It is also often the case that the homolog (be it orthologous or paralogous) of a protein in species B has a completely different function than its homolog in species A despite a high level of sequence similarity. This is usually very hard to determine in silico.
To find the functionally homologous protein(s), you would ideally need to identify the essential residues that allow your protein to perform its function. This could be done using something like PFam which will identify protein domains. You can then check whether the homologs you find also have this domain.
This is essentially what PSI-blast does. Although it doesn't take domains into account, each successive iteration is used to build a model of your proptein. The model is built under the reasonable assumption that highly conserved residues are important. So it will consider more diverged sequences as homologous if those residues are conserved.
If you know how your protein works and what residues are important, you can use that knowledge to refine the results of your PSI-blast. If you don't, you'll have to use only "good" hits to make the model. One way to do this, for well studied proteins, is to only add proteins that are already annotated as homologs of what you are looking for to build your model, then use that model to search in un-annotated species.
If you don't know, you could try looking for recognizable protein domains in your query protein (use PFam) and then use the HMM (hidden markov model) of the domain to identify important residues. For example, this is the HMM logo for the zf-A20 zinc finger domain:

The huge cysteine (C) residues are shown that size because they are very conserved across proteins carrying this domain and, presumably, are functionally important for the domain. So, if you pass your protein through PFam and identify domains, find the important residues and make sure all your homologs have those conserved. If using PSI-blast, only include sequences where those residues are conserved in the results you keep.
Finally, another useful tool that works in the same way is HMMER. This takes a protein alignment as input, like PSI-blast builds an HMM model from it and then can use this model to query a protein database for more hits. Methods like HMMER and PSI-blast are far better than simple sequence similarity approaches when looking for homologs.

Great tool HMMER, is as fast as BLAST and works pretty well

Dear terdon, your answer was the best answer someone could have given to this question and it was exactly the answer I was looking for. Thank you very much for your excellent answer

@Sara you're very welcome :)

It sounds like Smart BLAST might do what you want. Here's the description of it's goal:

SmartBLAST is a new and experimental NCBI tool that makes it easier to
  complete common sequence analysis tasks, such as finding a candidate
  protein name for a sequence, locating regions of high sequence
  conservation, or identifying regions covered by database sequences but
  missing from the query.
To do this, SmartBLAST performs the following tasks in much less time
  than it takes to run a typical BLASTp search:

a BLASTp comparison of the query with the closest matching sequences available;
a parallel BLASTp search to find the closest matches to high quality sequences from model organisms;
a multiple alignment between the query and five of the closest matching sequences (usually including two high quality sequences);
an analysis that produces a phylogenetic tree from the multiple sequence alignment.

[from NCBI Insights]

SmartBLAST is a very good idea as it gives so much more divergent organism, but is there any way to give its result to PSI-BLAST for more iterations?

@Sara just note the accession numbers of the sequences found by Smart BLAST and make sure to add them to your alignment in the psi-blast. They should be found as hits there too.


text:
How can we find the distance between all residues in a PDB file?
If we have a PDB structrure, how can we find residues physically interacting with each other in space? I know that we must find the distance between residues and if the distance is less than 5-6 Angstrom, we say that residues are physically interacting. But how can we find the distance between all residues and how can we finally determine the distances between all residues? Is there a software or webserver for that?

These interaction may be already annotated in the PDB, see Connectivity Section and Connectivity Annotation Section in the PDB spec.

One thing you're going to have to define for yourself is what you mean by "distance between residues". Do you mean closest atom pair? Closest heavyatom pair? Distance between C-alphas? The exact mechanism is probably going to be similar, regardless, but the details and the results will differ based on the definition.

I mean the distance between C-alphas

Is the Ring 2.0 webserver (http://protein.bio.unipd.it/ring/) good enough for this purpose?

If anyone is interested in speed benchmarks of distance calculations using various software, there are some at https://github.com/jgreener64/pdb-benchmarks. There are also example scripts for each software in that repository.

If you need to process multiple files, you could use Biopython to parse a PDB structure.
from Bio.PDB import PDBParser

# create parser
parser = PDBParser()

# read structure from file
structure = parser.get_structure('PHA-L', '1fat.pdb')

model = structure[0]
chain = model['A']

# this example uses only the first residue of a single chain.
# it is easy to extend this to multiple chains and residues.
for residue1 in chain:
    for residue2 in chain:
        if residue1 != residue2:
            # compute distance between CA atoms
            try:
                distance = residue1['CA'] - residue2['CA']
            except KeyError:
                ## no CA atom, e.g. for H_NAG
                continue
            if distance < 6:
                print(residue1, residue2, distance)
        # stop after first residue
        break

If you need to look at one structure, using a viewer perhaps would be easier. You could try PyMOL: (how to measure distance). There are other PDB viewers, some of which can work even through a browser.

Dear lakov, I installed Anaconda and Biopython packages and run the above code in PyCharm 5.0 and got the following error:

Traceback (most recent call last):
File "C:/Users/VN7-592/PycharmProjects/untitled/sara.py", line 7, in
structure = parser.get_structure('PHA-L', '1FAT.pdb')

File "C:\Users\VN7-592\Desktop\Anaconda\lib\site-packages\Bio\PDB\PDBParser.py", line 81, in get_structure
with as_handle(file, mode='rU') as handle:
File "C:\Users\VN7-592\Desktop\Anaconda\lib\contextlib.py", line 17, in enter
return self.gen.next()
File "C:\Users\VN7-592\Desktop\Anaconda\lib\site-packages\Bio\File.py", line 88, in as_handle
with open(handleish, mode, **kwargs) as fp:
IOError: [Errno 2] No such file or directory: '1FAT.pdb'

You don't have the 1FAT.pdb file in the working directory, or your operating system is case sensitive.

But this gives only intermolecular distances? BioPython ignores crystal symmetry? In general, to check for contacts between residues X and Y one would need to calculate distance between X and the nearest symmetrical image of Y.

Could you use CCP4's NCONT program? There's a GUI and a command line interface, whatever suits. You can specify which chains you want to target and interact with and set a cut off for distance. The bonus here is once you're in you have a nice suite of other structural tools to use.
If you're just doing it once, the GUI is friendly enough to work things out, if you're doing a batch then you can run it across several files via the command line.
Download

I have downloaded the program but I haven't been able to work with it yet. I'm trying to run it

@Sara what is it you're having a problem with?

It doesn't run with an error. The program says choose a project name and even when I have chosen a project name, again says choose a project name and stops running

As part of a project me and some teammates did a script that outputs visual maps of distances between residues. It uses Biopython.
The module contact_map.py does what you are looking for. As an example, if you want to find the residues whose CA are below 5 you can run the following command:
python3 contact_map.py pdb1cd8.ent -a CA -CA 5

This will produce three files:
distance_map_pdb1cd8_CA.png # Heatmap of the distance between the residues
contact_map_pdb1cd8_CA.png # Black/White heatmap: If it is at that min distance
contact_map.log  # The actions taken

If you don't have downloaded already the pdb structue you can use the main module:
python3 cozmic.py real 1cd8 -a CA -CA 5

How it calculates distances? Inside the asymmetric unit only, or as the shortest distance between the two residues in the crystal (or in biological assembly)?

@marcin It calculates the distance between the C alpha (or beta if you chose so) between each residue in the file.

Usually, the file == asymmetric unit. If you'd like to find all contacts in biological assembly or in crystal you'd need to take the symmetry into account.

Thanks for the tips, maybe you want to contribute to the modules.

I'm working now on a project that in the future may be competitive to Bio.PDB, but probably not this year yet :-)

I think Bio.PDB is being actively modified, but good luck with your work

You can use MDtraj. The package is easy to install using Anaconda.
You can get the interacting residues using the following snippet (taken from http://mdtraj.org/1.6.2/examples/native-contact.html)
heavy_pairs = np.array(
    [(i,j) for (i,j) in combinations(heavy, 2)
        if abs(native.topology.atom(i).residue.index - \
               native.topology.atom(j).residue.index) > 3])

# compute the distances between these pairs in the native state
heavy_pairs_distances = md.compute_distances(native[0], heavy_pairs)[0]
# and get the pairs s.t. the distance is less than NATIVE_CUTOFF
native_contacts = heavy_pairs[heavy_pairs_distances < NATIVE_CUTOFF]
print("Number of native contacts", len(native_contacts))

Distance between C-alphas can also be found using the BioStructures.jl package in Julia:
using BioStructures
struc = read("1AKE.pdb", PDB)
calphas = collectatoms(struc['A'], calphaselector)
dm = DistanceMap(calphas)

Then dm[5, 10] gives 12.19677604943208 etc.
Other variations also work, see the docs.


text:
Why is FPKM still used for gene expression studies?
Perhaps my understanding is wildly misguided, but I'm seriously confused about why people are still interested using FPKM values for cross sample gene expression analysis. My understanding is, that for some gene ${ i }$, its $FPKM$ is calculated as follows:
$$ FPKM_i \propto \frac{ \mathrm{fragments\ mapping \ to\ gene \ }i } { \sum_j \mathrm{ fragments\ mapping\ to \ gene\  } j } \bigg/
\mathrm{effective\_length}(i) $$
For comparing gene expression within a single sample, this doesn't seem too harmful to me. In fact, this seems like a good idea.
However, the use of $FPKM$ or even the famous $log_2(FPKM+1)$ in a similarity analysis seems to me plain wrong. True differences in RNA abundances of a couple of genes would significantly impact the FPKM values of the other genes leaving  perhaps confused about the true differences between the samples.
Am I missing something here? Is it ever meaningful to use FPKM values for analyzing across samples?

Is it ever meaningful to use FPKM values for analyzing across samples?

One should never use FPKMs for anything important. They can occasionally be useful for plotting, but even in that case one needs to construct the FPKMs from properly normalized data and not use the original definition.
It's becoming ever rarer to see people using FPKMs these days, thankfully.

Thanks, I also saw your response here https://www.biostars.org/p/317417/ as the answer to my followup.

So what is a proper alternative? TMM?

@LinuxBlanket TMM or another robust normalization method, yes.