How to find novel transcripts using GFFcompare? I am trying to find novel transcripts from an RNA-seq database. Based on the advice I got, it seemed that using Stringtie for transcript assembly is a good way to go, and it supports novel transcript discovery even with the reference GTF file provided (which apparently significantly improves the assembly process). So I tried to follow the protocol provided in the Nature Protocols paper that described the use of Stringtie; the paper also suggested using GFFcompare for comparing the assembled transcripts with the reference, to quantify and find the novel transcripts.

I followed the entire protocol, and this was the output of GFFcompare when comparing the reference GTF with that of one of the samples after mapping -

#= Summary for dataset: ./hisat2/ERR188044_chrX.gtf
#-----------------| Sensitivity | Precision  |
        Base level:    51.7     |    79.5    |
        Exon level:    46.7     |    85.2    |
      Intron level:    47.2     |    93.9    |
Intron chain level:    31.2     |    64.4    |
  Transcript level:    31.6     |    52.3    |
       Locus level:    36.6     |    50.1    |

     Matching intron chains:     580
       Matching transcripts:     664
              Matching loci:     397

          Missed exons:    4395/8804    ( 49.9%)
           Novel exons:     426/4874    (  8.7%)
        Missed introns:    3832/7946    ( 48.2%)
         Novel introns:      83/3992    (  2.1%)
           Missed loci:     610/1086    ( 56.2%)
            Novel loci:     273/795 ( 34.3%)

 Total union super-loci across all input datasets: 749
1270 out of 1270 consensus transcripts written in 188044_only.annotated.gtf (0 discarded as redundant)

In the paper however, the number of novel genes and transcripts is clearly mentioned -

So, I was wondering how they were able to calculate the number of novel transcripts and genes, from the data that GFFcompare provides (they even mention in the paper that they got it from GFFcompare). I see the numbers related to novel exons and novel introns, but how do I translate it to the number and identities of the novel transcripts/genes?

Cross-posted in Biostars. Let me know if I shouldn't do this.

The output of gffcompare includes several files per run (just like cuffcompare). Example for a run:

$ ls cuffcmp.* | sed 's/\t/\n/'
cuffcmp.combined.gtf
cuffcmp.loci
cuffcmp.output.gtf.refmap
cuffcmp.output.gtf.tmap
cuffcmp.stats
cuffcmp.tracking

From my experience, the easiest file to manipulate is the tmap file:

  Tab delimited file lists the most closely matching reference transcript for each query transcript. There is one row per query transcript, and the columns are as follows:

1     Reference gene name    The gene_name attribute of the reference GTF record for this transcript, if present. Otherwise gene_id is used.
2     Reference transcript id     The transcript_id attribute of the reference GTF record for this transcript
3     Class code  The type of relationship between the query transcripts in column 4 and the reference transcript (as described in the Class Codes section below)
4     Query gene id   The query (e.g., Stringtie) internal gene id
5     Query transcript id     The query internal transcript id
6     Number of exons     The number of exons in the query transcript
7     FPKM    The expression of this transcript expressed in FPKM
8     TPM     the estimated TPM for the transcript, if found in the query input file
9     Coverage    The estimated average depth of read coverage across the transcript.
10    Length  The length of the transcript
11    Major isoform ID    The query ID of the gene's major isoform
12    Reference match length  The length of the longest overlap with a reference, ‘-’ if there is no such exonic overlap

The important column to look at for transcript classification is actually the 3rd, "class code". The content is meaningful if you used option -r when running gffcompare.

In the tmap file you can access info on both transcript and gene ID directly, both from the assembler and the reference annotation.

Now what you need is:

n. assembled genes
n. novel genes
n. novel transcripts
n. transcripts matching annotation

At this point, it is essential to understand what the authors mean by "novel" (gene/transcript), does it mean "non previously annotated" (= class code u, for "unknown")? Also, what about "matching", does it mean "perfectly matching" (= class code =, for "complete, exact match of intron chain")?

Under these assumptions, you can do some easy stats with some bash help.

n. assembled genes:

$cat cuffcmp.output.gtf.tmap | sed "1d" | cut -f4 | sort | uniq | wc -l
n. novel genes:

$ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f4 | sort | uniq | wc -l
n. novel transcripts:

$ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f5 | sort | uniq | wc -l
n. transcripts matching annotation:

$ cat cuffcmp.output.gtf.tmap | awk '$3=="="{print $0}' | cut -f5 | sort | uniq | wc -l

I hope this helps. I think it's an easier approach than trying to "translate the numbers related to novel exons and novel introns to the number and identities of the novel transcripts/genes" from the stats file.

Although you have to think a bit what "class codes" to include in each count. Maybe my assumptions are too strict for your use case.


RNA seq fasta file annotation from alignment to reference matches I've got a fasta file with some RNA seq data and another csv file with the output from plast where I've aligned it to a reference using plastn. I'm struggling with figuring out a command to append my original RNA seq data fasta headers with the annotations.

I've looked at awk and sed but can't wrap my head around how to use them in this instance.

CSV output example:

Contig11    ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA   78.95   95  20  0   371 465 4627    4721    1e-14   82.4    75  625 1   0   15.2    0   5692    1   0   15.2    0

And example fasta contig:

>Contig11
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

How would I be able to parse through the csv to find a matching contig with the fasta file, and append the second csv column to the fasta header. Ideally to produce:

>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA....

Cheers for any help!

An awk approach (assuming your input "csv" file is space-separated):

$ awk '(NR==FNR){a[">"$1]=$2; next;} (/>/){$1=$1" "a[$1]}1;' foo.csv foo.fa
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

Explanation

NR==FNR : FNR is the line number of the current input file and NR is the line number of all input. When processing more than one file, the two will only be equal while reading the first file.
(NR==FNR){a[">"$1]=$2; next;} : if we're reading the first file, the csv, save the second field as the value in the associative array a, whose indices are a > followed by the 1st field. I use the > only to simplify matching it later in the fasta file. The next moves on to the next line so the rest of the command isn't run for the 1st file.
(/>/){$1=$1" "a[$1]}: Now we're in the 2nd file. If the line starts with a >, then set the 1st field to be itself, a space and then whatever was stored in the array a for this 1st field.
1; : this is awk shorthand for "print this line".

You could expand the above into:

awk '{
        if(NR==FNR){
            a[">"$1]=$2;
        }
        else{
            if(/>/){
                $1=$1" "a[$1]
            }
        }
        print;
      }' foo.csv foo.fa

If your files are small enough to be sorted easily, you can use the FastaToTbl and TblToFasta scripts I have posted before and then join:

$ join -o 1.1,2.2,1.2 <(FastaToTbl foo.fa | sort ) <(sort  foo.csv ) | TblToFasta
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
>Contig12 mamamama|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
AKHSALKJHDLKAJBDLKAJDBADADKJALDKJBAD


biopython-motifs , How can I create motif from sequences of different length? I want to create motif of different length kmers.
Following code works for kmers of same length but not different length.What to do ?
from Bio import motifs

sequences = SeqIO.parse("myfile.fasta","fasta")
    instances = []
    for data in sequences:
      instance = data.seq
      instances.append(instance)

motif = motifs.create(instances)

Motifs are not designed to be variable length in general since PSSMs/PVMs don't support insertions/deletions (indels). Instead, I suggest profile HMMs, such as those offered by HMMER3, which do support indels. You can parse HMMER3 output with Biopython via Bio.SearchIO.HmmerIO


gnu parallel for macs peak calling These are the list of files

Control_Input_sorted.bam_rem.bam
Control_H2BUb_sorted.bam_rem.bam
Control_IgG_sorted.bam_rem.bam
PTPN6_g2_6_H2Bub_sorted.bam_rem.bam
PTPN6_g2_6_Input_sorted.bam_rem.bam

So to get unique files i do this

 ls -1 *bam | sort | sed -r 's/_sorted.bam_rem.bam//g' | sort | uniq
Control_Input
Control_H2BUb
Control_IgG
PTPN6_g2_6_H2Bub
PTPN6_g2_6_Input

Now i have to run each sample against input and IgG for peak-calling.
Like  Control_H2BUb_sorted.bam_rem.bam against

Contol_Input_sorted.bam_rem.bam & Control_IgG_sorted.bam_rem.bam and

PTPN6_g2_6_H2Bub_sorted.bam_rem.bam against PTPN6_g2_6_Input_sorted.bam_rem.bam

So how do I parse and execute using gnu parallel.

I came across this tutorial where it does something like this

cat sample_names.txt | parallel --max-procs=12 'macs2 callpeak -t {}-A-NC.sorted.bam \
 -c {}-G-NC.sorted.bam -g hs -n {}-A-NC-sharp-model -q 0.01 --outdir {}-A-NC-sharp-model-peaks 2> {}-A-NC-sharp-model.stderr'

Now im not sure how to pass the file names as argument for IP and Input or IgG.

If i use logic a little bit i have to fix Control and PTPN6_g2_6 as constant naming!.

Any suggestion or help would be really appreciated.

not gnuparallel but nextflow. I don't know macs and what is 'first.bam' and 'second.bam' is not clear not me but I hope you'll get the idea.

The nextflow:

Channel.fromPath(params.bams).splitCsv(header: false,sep:'\t',strip:true).map{T->file(T[0])}.into{bams1;bams2}

boolean isControl(file) {
    return file.name.endsWith("_Input_sorted.bam_rem.bam") || file.name.endsWith("_IgG_sorted.bam_rem.bam ");
    }

bams1.filter{T->isControl(T)}.set{controls}
bams2.filter{T->!isControl(T)}.set{experiment}

process mac {
    maxForks 14
    tag "${bam1} vs ${bam2}"
    input:
        set bam1 , bam2  from experiment.combine(controls)
    output:
        file("output.peaks") into mac_peak
    script:
    """
    # check syntax please
    macs2 callpeak -t {bam1.toRealPath()}  -c ${bam2.toRealPath()} -g hs -n out-model -q 0.01 --outdir output.peaks
    """
    }

invoke:

find /dir1/dir2 -type f -name "*.bam" > input.list

nextflow run --bams input.list workflow.nf


Determine if a gene is mitochondrial or not I need to determine whether a gene is mitochondrial or not for C. Elegans automatically from its name by running a regular expression on a dozen of thousands of gene names. Currently I am having gene symbols like:

clec-190, clec-189, 21ur-8912, nduo-1, 21ur-6539, 21ur-13096, nduo-6, 21ur-14241, ndfl-4

I tried running grepl("^mt-", symbols) and it found nothing. I am sure that there are some mitochondrial genes just from their descriptions. I would assume that first I need to convert gene symbols to some other symbols to do that but I am not sure here.

Basically, I am trying to do it in a similar way to the tutorial of single cell RNAseq.

You can't determine if a gene is mitochondrial purely from its name. The simplest route would be to use the wormbase biomart and either download the gene names and their associated chromosomes or use the biomaRt package from bioconductor to do the same thing. You can upload a list of the genes you're interested in querying in either case (or just download all of them, and filter).

From your list, only nduo-1, nduo-6, and ndfl-4 are mitochondrial.


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

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.


How to interpret UniProt allele patterns? In creating a parser for the UniProt flat file uniprot_sprot.dat (directory) in the manual for the comments POLYMORPHISM there are listed patterns for alleles, e.g. for 1A11_HUMAN is

  The following alleles of A-11 are known:
  A*11:01 (A-11E), A*11:02 (A-11K), A*11:03, A*11:04, A*11:05 and A*11:07.

I tried searching for a description of the patterns using the words allele and/or polymorphism and found nothing.

How are the allele patterns to be interpreted?
Is there a search term to be used when looking for the patterns?
Also, is there an official repository of the alleles?

The answer does not have to explain the specifics of the example given.

Example of explanation of other allele pattern from here

This explanation of the patterns doesn't explain the part of the pattern in UniProt like (A-11E).

These are cytotoxic T-cell HLA alleles. HLA genotyping is very common and easy to do, so Genbank is the repository.

A*11:01 has very high frequency in Aborigine populations here . You can explore the population genetics and the past, present and current population genetics distributions per population of HLA at http://www.allelefrequencies.net/


How to know which ortholog database is the best for my work? I am working on four data sets of yeast genes and one data set of mouse gene. I need the list of human orthologs of those genes. There are a lot of databases available online which is in use for this purpose
https://questfororthologs.org/orthology_databases
I don't know which one to choose and why. Can anyone tell me on what basis should I choose a database over any other database? Thank you.

To recover orthologous sequence for model species
You should be able to find your orthologs using BioMart and your gene IDs.
There are numerous way to do it, directly on the ensembl website or through API. You can start from these links (quite old though) :

http://www.ensembl.info/2009/01/21/how-to-get-all-the-orthologous-genes-between-two-species/
https://www.youtube.com/watch?v=EwfOQsr2O-M #Get orthologs via BioMart
https://www.biostars.org/p/12245/ #Finding Human/Yeast Orthologs Using Embl?

Which database should I use
It is broad and depends on the nature of your work. In your case, working with "major model species", DIOPT (the first line of the table in Quest For Orthology) has the very feature you need on its front page for going from mouse to human. More eukaryotic/viral oriented database will be used for taxonomic identification, annotation inference, etc.
I hope it helps.


How to read the "SEQRES" section from a PDB file, using R I'm using R and its Rpdb package to work with PDB files. I use the read.pdb function but, nevertheless, I couldn't find the argument which will allow me to view the SEQRES section and the full amino-acid sequence of the protein.

For example, I want to find which amino-acid is in the 17th position of the sequence and I want to extract it from the 'SRQRES' section in the PDB file (the original length of the sequence was reduced for the sake of the example):

SEQRES   1 A  390  GLU PRO GLU TRP THR TYR PRO ARG LEU SER CYS GLN GLY
SEQRES   2 A  390  SER THR PHE GLN LYS ALA LEU LEU ILE SER PRO HIS ARG
SEQRES   3 A  390  PHE GLY GLU ALA ARG GLY ASN SER ALA PRO LEU ILE ILE
SEQRES   4 A  390  ARG GLU PRO PHE ILE ALA CYS GLY PRO LYS GLU CYS LYS
SEQRES   5 A  390  HIS PHE ALA LEU THR HIS TYR ALA ALA GLN PRO GLY GLY
SEQRES   6 A  390  TYR TYR ASN GLY THR ARG GLU ASP ARG ASN LYS LEU ARG
SEQRES   7 A  390  HIS LEU ILE SER VAL LYS LEU GLY LYS ILE PRO THR VAL
SEQRES   8 A  390  GLU ASN SER ILE PHE HIS MET ALA ALA TRP SER GLY SER
SEQRES   9 A  390  ALA CYS HIS ASP GLY ARG GLU TRP THR TYR ILE GLY VAL
SEQRES  10 A  390  ASP GLY PRO ASP SER ASN ALA LEU ILE LYS ILE LYS TYR

...

The original file can be found here 1a4g.pdb.

Any ideas on how to accomplish that?

I have not used Rpdb but I have used, and would recommend, bio3d.

To read the PDB file 1A4G:

library(bio3d)
pdb1a4g <- read.pdb("1a4g")

To get the SEQRES information, there are two options. Get it from the pdb object directly:

seq1a4g <- pdb1a4g$seqres

Or call pdbseq:

seq1a4g <- pdbseq(pdb1a4g)

Using the second method, the 17th position (note this is not the same as the 17th position of the original protein sequence):

seq1a4g[17]
 92
"Q"


Identifying relevant SNPs from a list I have a list of all SNPs (a list of RsIDs) that fall in the coding sequence of several thousand human genes, and I'm looking for a convenient method to come up with a list of known SNPs that alter an amino acid in those genes.

I'll be glad to hear your thoughts.

What you are looking for is SNP annotation. If you have the chromosome:position reference and alternate alleles for your SNPs of interest, it can be as simple as uploading them to the variant effect predictor.

http://grch37.ensembl.org/Homo_sapiens/Tools/VEP

This will give you the predicted protein change and novelty of the variant with respect to known variation in large scale sequencing projects like the 1000 genomes.

If you cannot afford to upload data on the web, there are downloadable command line software like SNPeff, Variant effect predictor etc but it would require some more bioinformatics knowledge of converting your data to a VCF format and obtaining a human reference genome.


How are Principal Component analyses and Admixture analyses from a genetic alignment different? How are Principal Component analyses and Admixture analyses from a genetic alignment different?

My understanding is that a PCA will take raw genetic differences across the entire alignment and plot them using dimensionality reduction techniques (single value decomposition for instance) whereas the Admixture analysis will plot likelihoods of sharing certain amounts of ancestry, based on allele frequency data. What I don't understand is how the two are fundamentally different. That is, isn't every variable position found in a PCA going to represent an allele used in Admixture?

Just to add to the great summary by Devon Ryan: admixture analysis tools are much more flexible than PCA (which is just a fixed mathematical operation), so they can be designed to incorporate LD patterns, phase information, different population evolution models etc.

Here's a paper from the creators of fineSTRUCTURE and Chromopainter admixture tools detailing the action of both - it seems they provide some comparisons with PCA as well:
http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002453


`TxDb.Hsapiens.UCSC.hg19.knownGene` and `EnsDb.Hsapiens.v75` can't be installed I think there is a problem with TxDb.Hsapiens.UCSC.hg19.knownGene and  EnsDb.Hsapiens.v75.  I'm struggling to install all of these packages as they pause at  byte-compile and prepare package for lazy loading . Can anyone try to re-install them in their local machine because am not the only one who is having the problem, please? They used to work fine until a week ago!
As usual, I install the packages by:
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

I also tried installing them individually:
BiocManager::install("EnsDb.Hsapiens.v75")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")

It pause at:
* installing *source* package ‘TxDb.Hsapiens.UCSC.hg19.knownGene’ ...
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading

sessionInfo( )
R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

loaded via a namespace (and not attached):
[1] BiocManager_1.30.10 compiler_4.0.3      tools_4.0.3

I've tried with and without RStudio.
I'd be grateful if someone can try to replicate the error, please?

I am able to successfully install these packages using:
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

This completes without any errors / warnings (R v4.0.3 / macOS 10.16), so I don't believe the packages themselves/Bioconductor are 'broken', but I suspect you are having compiler issues due to the Big Sur update.
These are the steps I used to solve the problem:

Reinstall xcode command line tools (even if it says "up to date")

sudo rm -rf /Library/Developer/CommandLineTools
sudo xcode-select --install

Install gcc & llvm via Homebrew (instructions for installing Homebrew) or, if you already have gcc/llvm installed via Homebrew, skip to the next step

# This can take several hours
brew install gcc
brew install llvm

Once you have gcc & llvm installed via Homebrew:

brew cleanup
brew update
brew upgrade
brew reinstall gcc
brew reinstall llvm

Link some headers into /usr/local/include

sudo ln -s /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/* /usr/local/include/

# You can safely ignore warnings like this:
#ln: /usr/local/include//tcl.h: File exists
#ln: /usr/local/include//tclDecls.h: File exists
#ln: /usr/local/include//tclPlatDecls.h: File exists
#ln: /usr/local/include//tclTomMath.h: File exists
#ln: /usr/local/include//tclTomMathDecls.h: File exists
#ln: /usr/local/include//tk.h: File exists
#ln: /usr/local/include//tkDecls.h: File exists
#ln: /usr/local/include//tkPlatDecls.h: File exists

Create a new ~/.R/Makevars file (if you already have a ~/.R/Makevars file, save it in a different directory (away from ~/.R/)) and include only these lines in the file:

FLIBS=-L/usr/local/gfortran/lib/gcc/x86_64-apple-darwin19/10.2.0 -L/usr/local/gfortran/lib -lgfortran -lquadmath -lm
CXX1X=/usr/local/gfortran/bin/g++
CXX98=/usr/local/gfortran/bin/g++
CXX11=/usr/local/gfortran/bin/g++
CXX14=/usr/local/gfortran/bin/g++
CXX17=/usr/local/gfortran/bin/g++

LLVM_LOC = /usr/local/opt/llvm
CC=/usr/local/gfortran/bin/gcc -fopenmp
CXX=/usr/local/gfortran/bin/g++ -fopenmp
CFLAGS=-g -O3 -Wall -pedantic -std=gnu99 -mtune=native -pipe
CXXFLAGS=-g -O3 -Wall -pedantic -std=c++11 -mtune=native -pipe
LDFLAGS=-L/usr/local/opt/gettext/lib -L$(LLVM_LOC)/lib -Wl,-rpath,$(LLVM_LOC)/lib
CPPFLAGS=-I/usr/local/opt/gettext/include -I$(LLVM_LOC)/include -I/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include

Compile a package from source in R/Rstudio

# To check whether openmp is enabled, compile data.table:
install.packages("data.table", type = "source")

If your package fails to compile, a couple of SO users have had to install a fresh gfortran (re: https://stackoverflow.com/a/65334247/12957340), which you can download from https://github.com/fxcoudert/gfortran-for-macOS/releases/tag/10.2-bigsur-intel

Once you can install packages from source without error/issue, update Bioconductor
BiocManager::install(version = "3.12")

Then, finally, you should be able to install your packages using
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))


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!

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.


How do I interpret my phylogenetic tree construction from Glucosidases in the Actinobacteria class? Dear knowledgeable bioinformaticians, I'm currently pursuing an undergraduate degree and as a result I do not have so much knowledge in this topic. I'm having some troubles understanding my phylogenetic three construction.

For the construction of a phylogenetic tree UGENE was used. 11 homologous proteins from different species in the same genus/phyla were aligned. For the tree construction PHYLYP Neighbour Joining (NJ) with the distance matrix model “Jones-Taylor-Thornton” was used. Two trees were constructed, one without bootstrapping and one with bootstrapping (1000 pseudoreplicates).

How do I interpret my values? Are they unusable since I have low bootstrap-values in a lot of my nodes? All guidance is appreciated, thank you.

Bootstrap values usually give you useful information, be they low or high. Usually, a "high" bootstrap value is something >= 70%. That means a clade is strongly supported; lower values are telling you the clade is not clearly separated from others. This could be for several reasons, e.g.:

1) A subpopulation is not really separated from others, e.g. the viral sequences in  town "A" and town "B" are completely admixed, meaning the population is panmictic and there is continuous gene flow between the two of them;

2) the populations diverged not long ago, it's too soon to have any phylogenetic signal, at least using the marker you are using (different genes evolve at different rates). Usually, housekeeping genes (which are usually very conserved) are used for populations that supposedly diverged long time ago.

The interpretation of your tree should take into account the biology and epidemiology of your target organism. If you really want to learn more, I suggest a lot of reading, and definitely a class dedicated to molecular phylogenetics.

This might be a good start:
https://www.amazon.com/Phylogenetic-Handbook-Practical-Approach-Hypothesis/dp/0521730716

Also, the NJ algorithm should allow for uneven branch lengths, maybe the tree you are showing is the consensus tree?

Hope it helps

Max


Get Salmon mapping/alignment summary With HISAT2, after the alignment of fastq files you get an alignment summary like this:
20000 reads; of these:
  20000 (100.00%) were unpaired; of these:
    1247 (6.24%) aligned 0 times
    18739 (93.69%) aligned exactly 1 time
    14 (0.07%) aligned >1 times
93.77% overall alignment rate

Is it possible to get this information with Salmon?

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


How to analyze co-occurrence of multiple SNPs? I am interested in 20 different SNPs that all are either As or Gs, and they all occur on the same chromosome. How can I assess the co-occurrence of these SNPs? In other words, I want to know, if SNP1 is "G", does that make SNP2 to be more likely to be a "G" as well?

I know this is broad and I'm not looking for a script or anything so specific, just a point in the right direction. I'm working with a non-model organism, so any organism-specific application won't be useful without being able to learn the methodology.

Googling thus far has only brought up GWAS and similar methods of correlating a genotype and a phenotype, but not within a genotype.

I found the LDheatmap package in R very helpful at calculating linkage disequilibrium. You can make plots such as this:


How to split FASTQ reads without re-running `fastq-dump`? I downloaded a few FASTQs from the SRA that were paired-end 76 bp.
When I look at the FASTQs, I get something like this:

@SRR7702334.1 1 length=152
AAGATNTTTTTCAGTATCAATCAGANGNNANANNCNNCTGCCCCAAATCCTGCAAAGAGCATGCAACGTCCTTAGAAGACAACGTCTTGCTCTGTCTCACCCAGGCTGGAGCATAGGGGCATGATCGATCATGACAAGGAGTTTACATTTTG

I forgot to add the --split-3 option to my fastq-dump command.
Is there an easy and fast way, without re-downloading the files, to split the reads in half and write them to two output FASTQs?

I tried the following using Biopython, and it'll work, but it's going to take a long time to finish.

from os.path import basename
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from gzip import open as gzopen
import argparse
from tqdm import tqdm

def split_read(read, pos=152 // 2):
    """
    Split read in half to produce 2 reads

    Parameters
    ----------
    read_id : Bio.SeqRecord.SeqRecord
        Read to be split
    """
    return {
        # R1/R2 as keys
        "R" + str(i + 1): SeqRecord(
            # first and second halves of sequence string
            seq=Seq(str(read.seq[(pos * i):(pos * (i + 1))]), IUPAC.ambiguous_dna),
            id=read.id,
            description="",
            letter_annotations={"phred_quality": read.letter_annotations["phred_quality"][(pos * i):(pos * (i + 1))]}
        ) for i in range(2)
    }

def split_fastq(fastq, prefix="split"):
    """
    Split reads from a FASTQ

    Parameters
    ----------
    fastq : str
        Path to input FASTQ file
    prefix : str
        Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)
    """
    # check for file format
    if fastq.endswith(".gz"):
        fq_handle = gzopen(fastq, "rt")
    else:
        fq_handle = open(fastq, "r")
    out_handles = {r: gzopen(prefix + "_" + r + ".fastq.gz", "wt") for r in ["R1", "R2"]}
    # load reads for random access
    records = SeqIO.parse(fq_handle, "fastq")
    for read in tqdm(records):
        # split read into each mate and fix header information
        mates = split_read(read)
        # write to each output file
        for m in mates:
            out_handles[m].write(mates[m].format("fastq"))
    fq_handle.close()
    for m in mates:
        out_handle[m].close()

if __name__ == "__main__":
    # parse command line arguments
    PARSER = argparse.ArgumentParser(description="Split a FASTQ into 2 FASTQs by splitting each read in half")
    PARSER.add_argument(
        "fq",
        type=str,
        help="Input FASTQ to split"
    )
    PARSER.add_argument(
        "-p",
        "--prefix",
        type=str,
        help="Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)",
        default="split"
    )
    ARGS = PARSER.parse_args()
    # run splitting function
    split_fastq(ARGS.fq, ARGS.prefix)

```

My comment above mentions the slow IO for the Python code, so I tried an alternative method that would cut down on the IO overhead. It's not super efficient, but I ended up using zcat | awk with a switch-case on the line number to stream the data.

# first pass to generate _R1 FASTQ
zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 1, 76); break; case 3: print "+"; break; case 0: print substr($1, 1, 76)} }' | gzip -nc > {output}_R1.fastq.gz

# second pass to generate _R2 FASTQ
zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 77, 76); break; case 3: print "+"; break; case 0: print substr($1, 77, 76)} }' | gzip -nc > {output}_R2.fastq.gz

All in all this took ~4 h to complete, which means an iteration speed of ~ 70 000 reads/s, or about 20x faster than the Biopython solution.
If I engineered a solution that would have only passed through the FASTQ once, it would've taken half the time, but in either case it was likely shorter than re-downloading the entire FASTQ.

But I will probably make use of fasterq-dump next time, and make sure I specify to split the read mates.


How can I use Nanopore reads to close gaps or resolve repeats in a short-read assembly? Low coverage MinION reads should be useful to close gaps and resolve repeats left by short-read assemblers. However, I haven't had any success with the software I know about. I'm aware of the following packages, either for scaffolding or closing gaps in short-read assemblies using long reads:

LINKS
OPERA-LG
npScarf
PBJelly
SSPACE-longread

I've tried npScarf and LINKS, and they ran successfully but didn't resolve the gaps in my assembly. I couldn't get PBJelly and OPERA-LG to run with the current versions of BLASR and samtools, and both packages seem not to be maintained. I have not tried SSPACE-longread because it's not open source.

What software can I use to fill gaps and resolve repeats in a short-read assembly with low-coverage Nanopore data?

More information

I'm finishing a mitochondrial genome. I have a 17 kb short read assembly with one gap. This was made from 2x150b paired-end reads from a TruSeq PCR-free library with a 400 b insert size. Molecular data suggests the genome may be around 32 kb.

After mapping short reads against the short read assembly, Pilon identified a tandem repeat of 156 bp in an AT-rich region, but wasn't able to close the gap.

I have 1.4 Gb of MinION rapid reads with an L50 of 8808 b. I mapped these reads against the short read assembly, and I can see reads that span the gap. I seem to have the information required to close the gap, but I don't know how to do it.

The organism has a nuclear genome > 600 Mb, so long read coverage is low. Despite this, I tried a de novo Canu assembly of the Nanopore reads, and I pulled out a 39 kb contig containing mitochondrial genes. Most genes on this contig are duplicated and fragmented, and Pilon wasn't able to improve it.

Thanks for reading!

You might want to look into Unicycler (manuscript with more information can be found here); even though it is supposed to be used with bacterial genomes only, it might work well with a small genome such as a mitochondrion's.

Beware that if you happen to have very long reads, you might end up with an assembly with multiple copies of the circular genome: you might want to look into circlator then.


How can I not show insertions in the Integrative Genome Viewer (IGV)? I am using IGV 2.5.2 and would not like to see the insertions  in my aligned reads. How can I do that? I have managed to remove mismatched bases because there is an option to do so when I right-click on the window ("Show mismatched bases"). Now I would like to do the same with these purple insertions.

3rd January 2020 update

Thanks to a member's response, I found some parameter in the preferences that might be helpful to hide the indels. However, I have tried different numbers from 0 to 1000 and the indels are still shown. What numbers should I put to the parameters "Show indel threshold (bases)" and should I click the box or not?

Preferences -> Alignments -> Hide indels < 'x' bases.

Credit to Pierre Lindenbaum


Filtering out reads from a reference (e.g. rRNA) using bowtie2 As I understand it, bowtie2 can easily be used to split reads into one of two groups:

reads for which both of a pair align well to a reference (using e.g. --al-conc-gz)
reads for which one or both of a pair do not align a reference (using e.g --un-conc-gz)

But I really want to split this second group into reads for which neither of a pair align to the reference.
bowtie2's report give a nice breakdown, but it's not clear to me how to extract specific members of these groups:
1000000 reads; of these:
  1000000 (100.00%) were paired; of these:
    958118 (95.81%) aligned concordantly 0 times
    1329 (0.13%) aligned concordantly exactly 1 time
    40553 (4.06%) aligned concordantly >1 times
    ----
    958118 pairs aligned concordantly 0 times; of these:
      179 (0.02%) aligned discordantly 1 time
    ----
    957939 pairs aligned 0 times concordantly or discordantly; of these:
      1915878 mates make up the pairs; of these:
        711518 (37.14%) aligned 0 times
        32134 (1.68%) aligned exactly 1 time
        1172226 (61.18%) aligned >1 times
64.42% overall alignment rate

In this case, I really want those that didn't align at all (i.e. the 35.58% not included in the 64.42% overall alignment rate which is the same as the 37.14% of the pairs that did not align concordantly).
Since in my case the reference is rRNA and many of those sequences are shorter than my reads, I'm not surprised that so there are so few concordant reads. I'm really interested in reads that aren't concordant and aren't discordant discordant either.
So, if I'm sticking to bowtie2, it looks like my only option is to take the resulting SAM/BAM file and filter what I want based on flags and then convert that back into FASTQ files.
Is there another option that I'm missing? It seems that the remaining output options relate to unpaired reads or SAM files.
I'm looking into other aligners now...
My goal is to go from FASTQ files to FASTQ files without extra steps, but I can go through a SAM/BAM file stage if I need to and even build a solution around that if needed.

If you want to remove rRNA, use tools that are made for this. We usually go for SortMeRNA, which can do exactly what you describe your goal is, i.e. go from FASTQ to FASTQ, splitting your reads into rRNA and non-rRNA.
In our experience, it is very accurate, if a bit slow at times.


A basic question about phylogenetic inference with Bayes-MCMC I am new to the world world of Bayesian-MCMC inference and have been working my way through various texts. I am having some difficulty understanding how MCMC works. I understand that it is a way of sampling from the posterior distribution. What is puzzling me, however, is how the probabilities of the posterior distribution have been calculated. That is, if the Markov Chain is sampling from the posterior distribution, that posterior distribution must have probability values on some axis. How were those calculated?

In short, the posterior is a probability value between 0 and 1 which is calculated from Bayes rule, which In phylogenetic usage is given as the product of the prior probability and the likelihood divided by a normalising constant. The posterior probability allows you to state the probability of some hypothesis given the data at hand (as opposed to the likelihood parlance, which is probability of data given hypothesis).

In a phylogenetic setting, the sampler works (roughly) like this: propose some value of a parameter (topology, branch lengths, etc.) and use the specified prior and calculated likelihood for that parameter value given the data to calculate the unnormalised posterior. Next, propose a new value, and repeat the same step. Compare the two values, and if the posterior probability of the second is greater, accept this parameter value and repeat the cycle again. In this way the Markov Chain moves through posterior space until it reaches a reasonably stable optimal approximation of the posterior distribution. Here, we want to get as many samples as possible to maximise the accuracy of our posterior. The maximum a posteriori (MAP) value of the parameter and its credible intervals, etc, will come from the distribution of these samples.


Use of heterozygous SNPs in cancer research: why? When reading about allelic fraction (AF) and SNPs in cancer research, they always mention the fact that they're using heterozygous SNPs (informative SNPs). Why is this? Why can't we use homozygous SNPs?

In this case, what is the reference base (from the reference genome)?

I’m no longer working in tumour sequencing so I’m by no means an expert. But in a nutshell, the reason is that, as indicated, homozygous SNPs aren’t informative: if your allelic fraction is 100%, meaning you only observe a single nucleotide at a given position, we don’t know whether we’re dealing with a tumour related mutation. In fact, we probably don’t because the sample likely contains healthy as well as tumour cells.

By contrast, tumour-specific variants always have an AF < 100%, in other words they exhibit some degree of heterozygosity.

In fact, in healthy tissues you’d expect to see one of three cases: your variant allele (= sequence base ≠ reference) has a frequency of either 0% or 100% (homozygous), or 50% (heterozygous), ignoring rare somatic mutations. With tumour samples, by contrast, the fraction will vary depending on the penetrance of the variant in the tumour, and on the fraction of healthy to tumour cells in the sample.


Are sequencing reads written differently when inserts and indexes are sequenced seperately? I was going over different sequencing library protocols and noticed that the adapter sequence can vary between protocols. In some adapter sequences the index sits between the primer binding site and the insert. In this case during sequencing the index will be sequenced and then the insert that comes right after it. In other cases, the adapter sequence is built so that the index is outside the primer binding site, i.e. index - primer binding site - insert (see two leftmost parts of the attached figure). Here the insert is sequenced using the Read1 primer and afterwards the index is sequenced in a separate reaction using the 17 index primer.

In this case are the index reads just added on to the ends of the insert reads in the FastQ files? Is there a way to tell from the FastQ files where the index was positioned within the adapter?

The order of sequencing in a dual index flow cell is:

Read 1
Index 1 (I7)
Index 2 (I5)
Read 2

Between each step there is a switch of the primer used (Read1-primer -> Index1-primer (I7) -> Index2-primer (I5) -> Read2-primer)
There are differences in how the I5 is dealt with between dual- and single-read flow cells as explained in the Illumina Indexed Sequencing Overview. In this document there is a lot more detailed explanation of the sequencing workflow that may be helpful to you.
But as explained by the other answers, none of this will matter in most cases. Even if you decide to output the Index sequences they will be exported as separate FASTQ files.


How do population genetics people define a population? How do population genetics people define a population?

Do they define it as a layman will do? say Africans, Americans, or so? Or is there a more scientific way of doing so? For example, I think defining one's population as one's allele frequencies makes more sense, but has this point been discussed previously?

A population is defined as a group of people (or organisms of the same species) living together in a geographical region. However, for human genetic studies, it is wise to use ethnicity as well (so group of people with the same ethnic background). See here how they define population in the 1000 genomes project.


Correlation network representation I have troubles understanding the figure 10 of this paper:

  Similarity network of participating methods for BPO. Similarities are computed as Pearson’s correlation coefficient between methods, with a 0.75 cutoff for illustration purposes. A unique color is assigned to all methods submitted under the same principal investigator. Not evaluated (organizers’) methods are shown in triangles, while benchmark methods (Naïve and BLAST) are shown in squares. The top-ten methods are highlighted with enlarged nodes and circled in red. The edge width indicates the strength of similarity. Nodes are labeled with the name of the methods followed by “-team(model)” if multiple teams/models were submitted

This study calculated the pairwise Pearson correlation between methods on a common set of gene-concept pairs and then visualized these similarities as networks. How is the place of the methods determined in this visualization?

For example, the relationship between Pannzer (1) and Pannzer (2) is greater than the relationship between Pannzer (1) and Pannzer (3), but the Panzer (1) and Panzer (3) are closer to each other than Panzer (1) and Panzer (2). Is this proximity positioning random?  or is there a certain rule?

When plotting networks where only the weight of an edge is known like this one. The relative position between elements of the networks is irrelevant. The only informative part of the graph is (usually) the size of the edges, where is the thicker means more weight. This isn't a PCA, so the distance between different nodes is not relevant.

For those reasons there are several methods to build these networks, some of them give more emphasis to represent better the degree of the network or other characteristics of the networks.

I consider this representations bad for informing about the networks, because the cut off is arbitrary and the distance between nodes is not informative. One could first do a PCA or MDS and then add the edges between the nodes, this way the relative distance between the nodes would be meaningful as well as the edges. Or one could also use a heatmap to represent all the correlations between the nodes.


Pathway level analysis of single-cell gene expression I'm looking for single-cell specific methods to construct (using gene expression data) new features that express pathway "level" or "activity", and then use these for clustering cells.

One example for bulk RNA-seq is PLAGE, implemented in the GSVA R package.

I have found only one method for single-cell RNA-seq that comes close: pagoda in the scde R package. This performs PCAs for each pathway, using its genes.

A very similar question is Classification (supervised learning) of expression data on pathway level, on bulk/microarray data.

Package Hipathia allows you to compute a value of "activity" for each pathway and each cell, so that the matrix of gene expression is transformed into a matrix of pathway activity. You can use this matrix for further analysis, cell differentiation and clustering. You can also compute functional activity matrices in the same way.


What is the correct way to map Hi-C data with bwa mem? Library Prep

I have a Hi-C library prepped using an enzyme that cuts at GATC, so it leaves GATCGATC as the junction sequences. This library was sequenced on a 2x150 PE Illumina run.

Data Pre-processing

The reads were adapter and quality trimming with trimmomatic paired-end mode, keeping reads that overlap.
Finally, the forward and reverse reads were trimmed such that all bases 3' of the junction sequence (GATCGATC). For example, the read:

5' AAAATT**GATCGATC**AATTATGGATAA 3' becomes: 5' AAAATT**GATC 3'

In this process, sequences that don't contain the junction sequences were left alone. Trimming junction sequences in this way makes it such that each read pair maps to only two places in the genome. No three-way or four-way chimeric mappings are possible since everything 3' of the junction sequence was removed.

My Question

What is the correct way to now map these trimmed Hi-C reads to a genome using bwa mem?

Progress so far/things I have tried

In the past I have just tried mapping using bwa mem with no special parameters using this command: bwa mem ref.fa hi-c_f.fq.gz hi-c_r.fq.gz. One concern is that the computations that bwa mem performs on insert sizes in this mode will cause some mappings to be faulty or that read pairs with chromosome-length information will be removed from the alignment because they fall on the tail of the insert size distribution (non-Gaussian, but log decay).
Another option to consider is mapping each file, hi-c_f.fq.gz and hi-c_r.fq.gz, independently then adding mate information with samtools fixmate. However, I do not know how to do this.
A last option that I have considered is to use bwa mem in paired-end mode, but to use a special flag. In the past I believe that -5 was the correct flag, but in this post Heng Li discusses that the -5 flag behavior has changed, as well as a new -q flag. In this post Heng Li discusses that he thinks bwa mem is still a good option for mapping Hi-C data but does not discuss the proper flags.

Limitations

The answer should provide the correct way to map Hi-C using bwa mem, but other software suggestions will also be considered.

See here, in particular slide #10 from the tutorial:

bwa mem -5SPM ref.fa read1.fq read2.fq > out.sam

Here, -SP disables pairing. The Aiden lab and this paper also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.


Tools for creating gene coexpression networks What are some good tools (in Python if possible) for creating gene co-expression networks from RNA seq data?

After looking around, it seems like the Camoco library is suitable.


Error using htseq-count: Could not retrieve index file I have a total RNAseq dataset that I aligned using STAR producing BAM files (sorted by coordinates). I am now trying to get counts for the lncRNA sequences using htseq-count, with the command:
htseq-count -m intersection-nonempty -f bam -i gene_id -r pos -s reverse ./OVCA429_Aligned.sortedByCoord.out.bam ./GRCh38_latest_genomic.gtf > ./OVCA429.counts.txt

This is the output:
[E::idx_find_and_load] Could not retrieve index file for './OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam'
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
...(truncated output here)...
3775001 GFF lines processed.
[E::idx_find_and_load] Could not retrieve index file for './OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam'
100000 alignment records processed.
200000 alignment records processed.
300000 alignment records processed.
...(truncated output here)...

Output text file:
__no_feature    20368875
__ambiguous 0
__too_low_aQual 0
__not_aligned   0
__alignment_not_unique  2634439

I've read through the htseq-count literature but it never asks to include an index file, therefore i'm not sure how to deal with this problem. Any help/ suggestions would be welcome! Thank you, in advance.

The problem seems to be missing index files of the bam files (therefore the error message that the index file is missing). You should be able to generate those with samtools, e.g.
samtools index ./OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam


Can't install newest Blast from Conda I have run into an odd situation, where Anaconda repo page clearly states that the newest version of Blast in Bioconda channel is 2.7.1, however the version Conda wants to install (command conda install -c bioconda blast) is 2.6.0.

Even specifying the version with command conda install -c bioconda blast=2.7.1 does not work and Conda returns a PackagesNotFoundError.

By observing the webpage mentioned above, one can see that there is a "v" before the version number, however adding this letter to the command with specified version did not help.

To make things even more odd, conda search blast returns a list of Blast versions available and there is a 2.7.1 version available in Bioconda channel.

I have added the Bioconda channel into my Conda.

My question is: How can I install the 2.7.1 (by the time of writing this question the newest version available in Anaconda repo) through Conda?

Edit: I have Ubuntu Server 18.04 LTS and Conda 4.5.11.

I was able to reproduce your problem, even in a blank environment.

$ conda create -n test
$ conda activate test
$ conda install bioconda::blast=2.7.1
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - bioconda::blast=2.7.1
  - boost=1.64
  - bioconda::blast=2.7.1
  - gnutls
  - bioconda::blast=2.7.1
  - nettle[version='3.3|3.3.*']
  - bioconda::blast=2.7.1

Note that boost, gnutls, and nettle are listed as requirements here.

If you do a search for those packages, you'll see that boost is available in a default channel, whereas gnutls and nettle are only in the conda-forge channel, which is not a default.
Adding those packages to the install command resolved this issue, for me.

$ conda install bioconda::blast=2.7.1 conda-forge::gnutls conda-forge::nettle --dry-run
Solving environment: done

## Package Plan ##

environment location: /Users/hawleyj/anaconda3/envs/test

added / updated specs:
    - bioconda::blast=2.7.1
    - conda-forge::gnutls
    - conda-forge::nettle
...

As mentioned in the comment below, to avoid this issue entirely, before using Bioconda packages make sure you run the proper channel configuration first.

conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge


How to Convert From Gene Location to SNP Location I am new to this. I have an SNP array composed of digits [0,1,2] for chromosome 21. How do I find the start:end locations for the region that is "21q22.11" step by step?

You first need to get the coordinates of the "cytoband" 21q22.11, which I would search for using the UCSC table browser. You need the following settings:
Clade: Mammal
genome: Human
assembly:
group: Mapping and Sequencing
track: Chromosome Band
table: cytoBand
region: position chr21
output format: BED

(you can leave other options as default)

That will give you a bed file of cytobands on chr21. You'll need to search (e.g. grep) for your cytoband: q22.11.


differential analysis of chip-seq data I have several sets of chip-seq data. I called the peaks using Macs2. I am pretty new to the field and I will appreciate any help. I wanted to annotate the peaks and see which peaks are shared between them and which ones are different. What is the best way to do this? which package should I use?

Many thanks

For sharing, you can test the proportion of overlapping peaks using bedtools intersect or the find.overlap function in GenomicRanges packages in R.
What makes more sense is perhaps to merge all the peaks together using bedtools merge, and then count the number of reads/fragments in each peak for each sample using featureCounts. The output of this is a matrix, which is very similar to gene expression results. Then you can do some differential tests on this data.


Why do reverse radtags have different start points in radtag sequencing? Why is it that when doing paired end Radseq that R1 and R2 have different start points? For instance, while doing Radseq you do a restriction digestion of your DNA using enzyme X- shouldn't the sequencing reaction take place from either end of where the enzyme cut? If it doesn't- but happens randomly where do your multiplex sequences go?

For instance- if you are sequencing the sequence on an illumina machine with the "<>" representing the cut site-

Where does the sequencing reaction start from in each direction?

Where does the multiplexing sequence go? - does it go on each?

5-GGCTAGC<CCTGCA*GG>GCGCTT-3
 3-CCGATCG<GG*ACGTCC>CGCGAA-5

Assume that after digestion you're left with the following sequence:

5' AGCTAGCTAGCTAGCTGACGTCGTAC 3'
3' GTACGACGTCAGCTAGCTAGCTAGCT 5'

Then one read would start with AGCTAG and the other would start with TCGATC.

The barcodes (what are used for multiplexing) are outside of this and part of the illumina adapters. They're sequences in completely different reads. This is nicely illustrated below:

Index 1 and 2 don't actually need to be present. Even if they are, whether the sequencer actually sequences them is a user setting.

For those curious, the order of the reads as they're sequenced are:

read #1
Index 1 (optional)
Index 2 (optional)
read #2 (optional)

If you ever get more than 2 reads from a sequencing provider, then it's likely that they included the index sequence (I have no clue why, but some sequencing centers like to do this...even though it confuses the heck out of their users).


Mini-bulk RNA-seq and scRNA-seq analysis differences? We have RNA-seq data from libraries prepared using the Smartseq2 single-cell protocol on 500 cells (mini-bulk) / library. The complexity is much better than with single-cells (~14k genes for 1.5M reads).

6 cell types / biological replicates were collected with flow cytometry, and there are 4 control and 5 diseased biological replicates (patients).

We are mostly interested in the differences between control and disease either overall or within each cell type.

Can single-cell tools (e.g. Seurat) be used on this dataset? It is neither bulk nor single-cell data. How would you approach this?

Single-cell tools were developed to address some specific factors unique to scRNA-seq, such as:

a lot of dropouts, so hundreds or thousands of detected genes, not tens of thousands (for mammalian genomes)
hundreds or thousands of replicates
not necessarily well-defined groups

Thus, I would say your experiment is a lot closer to bulk RNA-seq than scRNA-seq.


remove sequences from fasta file matching a string in the header I have a file with 16S sequences. some headers contain species information. For my purposes I would like to exclude a number of species from the file, therefore I would like to do a pattern matching within a 16S header (that contains other characters beyond the species name) to then remove that header and the following sequence.

So, I have:

>NCBI_ID1;high_rank_taxonomy;...;species_name1
ATCCTT
>NCBI_ID2;high_rank_taxonomy;...;species_name2
ATTTTCCC
>NCBI_ID3;high_rank_taxonomy;...;species_name3
TTTCCCCAAAA
>NCBI_ID4;high_rank_taxonomy;...;species_name4
CCCTTTAAA

I would like to exclude, say, species 1,2 and 4 to retain:

>NCBI_ID3;high_rank_taxonomy;...;species_name3
TTTCCCCAAAA

I found this solution in an adjacent question:

awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa

where I have my species names in the file ids.txt and sequences in file seq.fa

but that would work if the whole of my header was just the species name (contents of each line in ids.txt), which is not the case.
I thought I could be smart and use *species_name1* instead of species_name1 as contents of ids.txt but that apparently doesn't work.
Could someone help me with modifying the current solution to accommodate the existence of additinal strings in the headers?

Many thanks!
A

you may use simple R commands, using the 'seqinr' package

#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("filename.fa")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), ";"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c("species_name1", "species_name2", "species_name4")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fa")


GenomicRanges: get nearest neighbor distances for random genes using a for loop I want to compute the nearest neighbour distances using a for loop. I can do this with a random subset by doing this:

df <-  read.table("tmp.txt", header=F)
colnames(df) <- c("gene", "seqnames", "start", "stop", "LOCID")
gr <- makeGRangesFromDataFrame(df, keep.extra.columns=TRUE)

# genes of interest in a list
gene.list <- read.table("careplastic_diffs.txt")

# filters the GenomicRanges object with genes of interest
sub.gr <- gr[(elementMetadata(gr)[, "gene"] %in% gene.list$V1)]

# compute distances
dists <- distanceToNearest(x = sub.gr)

mcols(dists)[,1] # contains a list of distances

I now want to essentially do this 5 times, using the same logic as above:

## create a rand dataframe, each column is a random list of 50 genes
rand.df <- data.frame(replicate(5, sample(gene.list$V1, 50)))

# loop through columns and gets distances and prints to a file
for(i in 1:ncol(rand.df)){
  rand.gr[i] <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists[i] <- distanceToNearest(x = rand.gr[i])
  output <- print(mcols(dists[i][,1]))
}

But I can't seem to make this work.

  Error: subscript contains out-of-bounds indices In addition: Warning
  message: In recycleSingleBracketReplacementValue(value, x, nsbs) :
  number of values supplied is not a sub-multiple of the number of
  values to be replaced

It's probably something really silly but any help would be appreciated ...

Is it because rand.gr has not been initialised, so when you ask to set rand.gr[i] it can't because there isn't one yet?

If you set i to a value (e.g. 1), and run the lines inside your loop, which of the lines produces the error?

If you don't need to retain all the random genomic ranges at the end, you could just leave out some of the [i]s, something like this:

for(i in 1:ncol(rand.df)){
  rand.gr <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists <- distanceToNearest(x = rand.gr)
  output <- print(mcols(dists[,1]))
}


Ranking metric in GSEA I am using fgsea in r to calculate and plot a bunch of GSEA graphs.

My question is whether it is acceptable to use the Wald statistic from DeSeq2 to rank the gene list?

I have seen in the GSEA application that the signal to noise ratio is used but is the Wald statistic sufficient?

If not then is there an easyish way to calculate the signal to noise ratio in R given a normalised count matrix?

Thanks in advance.

If using DESeq2 with GSEA, I'd recommend ranking by shrunk log2FC values.

It'd also be worth considering ranking positive and negative associations separately, because the standard GSEA algorithm doesn't cross at the zero point when associations change from positive to negative.

You shouldn't be using p-values to rank anything. A p-value gives information about the reliability of a result, but not its importance. See point #5 of this article about p-values:

https://amstat.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#_i31


How to convert a SDF molecular structure file to a pandas dataframe with Python? What is the best python tool to convert a SDF file to a structured data frame (pandas.DataFrame). Similar to Molconvert from ChemAxon? SDF can have different datafields for each compound which I would like to have in a dataframe. The 2D or 3D structures can be disregarded.

It is easily done with RDKit:
from rdkit.Chem import PandasTools

fn = '/mypath/dataset.sdf'

df = PandasTools.LoadSDF(fn, embedProps=True, molColName=None, smilesName='smiles')


Which quality score encoding does PacBio use? Do you know which quality score encoding PacBio uses now? I know some of their file formats have changed in the past year or two, but I haven't found much on their quality score encoding.

The most recent answer I found is from 2012, where one user states it's "phred style", but doesn't state whether it's 33 or 64: DeNovo assembly using pacBio data

PacBio does use PHRED 33, but it turns out the question may be irrelevant for the newer PacBio Sequel Sequencer, because it reports all base qualities as PHRED 0 (ASCII !). The RS-II reports PHRED 33 quality scores.

The Sequel provides data in the form of subreads, which are the circular consensus sequences (CCS) from a single zero-mode waveguide (ZMW). Essentially, a single circularized DNA molecule enters a ZMW well, and is sequenced as many times as possible. (Sometimes multiple molecules get into a single well, but that's not the normal case.) The number of times this single molecule can be sequenced depends on several variables, including molecule length. The consensus sequence for this single molecule is then reported out as a subread.

As I started playing with the quality scores in my data, I realized the quality score for every base is !, which equates to a PHRED quality of 0. I contacted someone at PacBio who informed me that they learned the base quality values were not any more useful than the bases themselves for generating a high-quality circular consensus sequence. Thus, rather than spend the computational effort to calculate the scores, they skipped it, and reported ! for every base. They also report an overall read quality of 0.8, but this is simply a placeholder. They are working to make this information more readily available.


Subsetting of object existing of two samples I have a Seurat object existing of an aggregate of two samples namely; RD1 and RD2. I am trying to make a subset of each sample. But it's not really working. I manage to make a subset of the cells but I don't manage to make a subset of the whole object. Isn't there a way to do this, based on the subset from cell names a subset to make from the whole object.

seuratObj$orignal_1 <- seuratObj@meta.data$orig.ident
seuratObj$orignal_2 <- seuratObj@meta.data$orig.ident
seuratObj$orignal_1 <- gsub("RD1", "TRUE",seuratObj$orignal_1 )
seuratObj$orignal_1  <- gsub("RD2", "False",seuratObj$orignal_1  )
seuratObj$orignal_2  <- gsub("RD1k", "TRUE",seuratObj$orignal_2  )
seuratObj$orignal_2  <- gsub("RD2", "False",seuratObj$orignal_2  )

subset_DR1 <- data.frame(orignal_1 = seuratObj$orignal_1)
subset_DR1 <- subset(subset_DR1, orignal_1 == TRUE,select = c("orignal_1"))

subset_DR2 <- data.frame(orignal_2 = seuratObj$orignal_2)
subset_DR2 <- subset(subset_DR2 , orignal_2 == TRUE,select = c("orignal_2"))

SubsetData(seuratObj, cells = rownames(subset_DR1), subset.name = seuratObj$orignal_1)

I just started to work with seurat packages. I think I am missing something.

I think your issue is simply not assigning the newly-subset object in to a new object.

Assuming rownames(subset_DR1) is indeed a list of cell names

seuratObj_subset_dr1 <- SubsetData(seuratObj, cells=rownames(subset_DR1))

will work.


Better aligner than bowtie2? Bowtie2 is probably the most widely used aligner because of it's speed. Burrow-wheeler (BW) algorithms (including bwa) tend to be faster. However, they have limitations when it comes to aligning very short reads (e.g. gRNA). Also, setting maximum number of mismatches allowed is complicated by the seed length, overlaps and other parameters.

I wonder if there is any better multi-purpose aligner out there. May be with algorithm other that BW. One which allows special cases e.g. allowing shord reads and high number of mismatches.

Bowtie2 is no longer the fastest aligner. Salmon and Kallisto are much faster, but have been designed to optimise RNASeq mapping. Their speed is gained from avoiding a strict base-to-base alignment, but they can output mostly-aligned reads (i.e. position-only, without local alignment) as pseudo-alignments. See here for more details.
Both Kallisto and Salmon can do additional bootstrapping that is interpreted by sleuth for improved performance in isoform detection. They can also output counts that are equivalent to read-level counts from other programs, which can then be used by other downstream gene-based differential expression analysis software (e.g. DESeq2). Salmon has additional options that can  correct mappings for sequence-level and GC bias.
HISAT2 is from the same group as Bowtie2, and does the same sort of stuff, but with a few optimisations added on top. In particular, it's much better at working out split reads from RNASeq runs. Like Bowtie2, it will do local alignment of reads.
For quick alignment of long reads, minimap2 works well. For high-accuracy alignment (but comparatively slower), LAST works well.
Most bioinformaticians seem to prefer STAR for things that Bowtie2 was previously used for. I'm not yet convinced it's a better alternative, and currently prefer HISAT2 for high accuracy short-read alignment.


Calculating Unifrac Distance and missing data I am trying to calculate unifrac distance for my OTU table, which has about 3000 taxa, and 300 samples. My phylogenetic tree only considers the most important 500 taxa. My questions are:

Do I ignore the 2500 taxa that were not considered in the phylogenetic tree?
Is their a package in R that deals with this kind of missing data?

To calculate the unifrac distance you need the phylogenetic tree. If you only have it for 500 taxa then you can only calculate it for those taxa.

It is impossible to impute the phylogenetic tree without knowing the sequence on which the taxa is based. But if you know it then you can calculate the phylogenetic tree for the 3k taxa you have.


How to download SNP data from specific regions, population and positions? I am totally new in Bioinformatics and I would like to apply my knowledge in feature selection on the tag SNP problems.
To do that, I've read a lot of papers and books in order to understand the main concepts and how the algorithms to this task work.
However, is been very hard for me to understand the organization of the huge amount of genomic data available over the internet.. HapMap, 1000 Genomes, dbSNP and so on... (!) I downloaded GB of data from different sources, but I really do not know what to do with then.
--
To be more specific in my question:
First I would like to run a Genetic Algorithm to find the tag SNP based on a simple linkage disequilibrium coverage measure. The authors of this GA used the following data

Three published SNP data sets were downloaded from HapMap (http://www.hapmap.org/) (NOT AVAILABLE ANYMORE) for evaluating the prediction accuracies achievable by our method:
• ENr112. 343 SNPs genotyped in population ASW (African ancestry in Southwest USA) on chromosome 2 from position 51512208 to 52012208. The genotypes of 83 individuals were used.
• ENr113. 380 SNPs genotyped in population ASW on chromosome 4 from position 118466103 to 118966103. The genotypes of 83 individuals were used.
• ENm013. 683 SNPs genotyped in population ASW on chromosome 7 from position 89621624 to 90736048. The genotypes of 83 individuals were used.

Besides that, the authors stated that

We represent a genotype g_i of length L by a sequence over {0, 1, 2}^L (since we assume only bi-allelic SNPs). 0 and 1 stand for the minor and major homozygous types respectively, meaning that the two alleles at that SNP locus are either 0 or 1; the major and minor alleles are determined by the frequency at which each allele appears in the population. A value of 2 denotes the heterozygous type for an SNP having two different alleles.

And also

To obtain phasing of genotypes, we used the GEVALT algorithm.

I believe they obtain the aforementioned data in genotype format (something with 0, 1 and 2 values only) and then phasing it to apply in the problem. The image below is a example of this process.

So

How and where can I download these data?
I know HapMap is not available due to a security issue but there are the FTP, but the FTP folders are really complex to understand. Can someone guide me to download data from there?

I hope I made myself clear.

The HapMap website was suspended due to a security issue. Archived HapMap data is available via FTP, but it's a better idea to download up-to-date International Genome data from the IGSR website (which includes the HapMap individuals as a subset).

The format you use will be dependent on the data analysis program you're using. In most cases involving SNP analysis, the VCF format is probably what would be used.


public multi-modal single-cell data There is a scRNAseq Bioconductor package with a few different example scRNA-seq datasets. Are there any R packages that offer multiple modalities of single-cell data? For example, hashtags or ADTs or ATAC-seq. I am trying to find something that can be used for testing or for a quick tutorial. Preferably something that is on CRAN or Bioconductor.

These tutorials on Seurat multimodal data and the wrapper Seurat data are easy ways to start.

The wrapper has some cite-seq data preinstalled making it easy to work with benchmarked data sets

If you are using Seurat u can just as well have a look at Signac to start working with some sc-atac seq


How to obtain desired output? I am working on a project using the following command within nano:
from Bio import SeqIO
import sys
import re

    fasta_file = (sys.argv[1])
    for myfile in SeqIO.parse(fasta_file, "fasta"):
      if len(myfile) > 250:
       gene_id = myfile.id
       mylist = re.match(r"H149xcV_\w+_\w+_\w+", gene_id)
       print (">"+mylist.group(0))   # edited from list to mylist

and its providing with the following outout:
>H149xcV_Fge342_r3_h2_d1
>H149xcV_bTr423_r3_h2_d1
>H149xcV_kN893_r3_h2_d1
>H149xcV_DNp021_r3_h2_d1
>H149xcV_JEP3324_r3_h2_d1
>H149xcV_SRt424234_r3_h2_d1

How can I change my command so that it provides me with the desired format and only the UNIQUE gene id's:
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2

its easy replace match with sub but please stop using list as a variable ... mylist is fine.
This might work
   mylist = re.sub(r'H149xcV_\w+_\w+_\w+', gene_id)

otherwise
myregex = re.compile('_\w+\s+.*')
fastaid = myregex.sub('', myfile)

OR from @MaximilianPress
myregex2 = re.compile('_\w+\n') # or myregex2 = re.compile('_\w+$')
fastaid2 = myregex.sub('\n', myfile) # or fasaid2 = myregex.sub('', myfile)

The above will work ... As with all my code I never proof it ..


Find number of possibilities We have a column which has following format : (1, 1, 1, 0, 1).

Here is the definition of what means two columns are compatible.

  Using the notation Oi to denote the collection of rows possessing a 1
  in their i-th element, we conclude that one of three possibilities
  must be true: Oi ⊆ Oj , Oj ⊆ Oi , or Oi and Oj are disjoint (the first
  two cases include the possibility that Oi = Oj). If columns i and j
  satisfy this condition, then we call them compatible

The question is : How many possible columns of length 5 are compatible with the column (1, 1, 1, 0, 1)?

You can read more here where problem is described in page 14-16.

Thanks in advance !

Let $n$ be the number of rows (in your case $n=5$) and $|O_j|$ the cardinality of $O_i$.

We can calculate the number of compatible columns as sum of the number of $O_i \subseteq O_j$, $O_j \subseteq O_i$ and $O_i, O_j$ disjoint.

number of $O_i \subseteq O_j$:
Here we need to find the number of supersets of our $O_i$. That means we need to find the number of sets which contain the indices in $O_i$ and for each index not in $O_i$, we have two choices: It is contained, or it is not contained in $O_j$. Therefore the number of supersets is $2^{n-|O_i|}$, which is 2 to the power of elements not in our set.

number of $O_i, O_j$ disjoint:
Here all elements of $O_i$ are fixed not to be in $O_j$ and for each index not in $O_i$, it can again be in $O_j$ or not. That means that the number of possible sets is again $2^{n-|O_i|}$.

number of $O_j \subseteq O_i$:
Here we need to find the number of subsets of $O_i$ and this is the number of elements in its power set, which is $2^{|O_i|}$.

So finally, this means that the number of compatible columns can be calculated as: $2^{n-|O_i|} + 2^{n-|O_i|} + 2^{|O_i|} - 2$.

Why -2? We have counted 2 sets twice:

The empty set, which is a subset of $O_i$ and also counted in the number of disjoint sets
and $O_i$, which we counted in $O_i \subseteq O_j$ and in $O_j \subseteq O_i$

This means for your examle, we have $16 + 2 + 2 - 2 = 18$ compatible sets.

Note: The all zero vector (which corresponds to $O_i=\{\}$) is a special case for which the formula above does not work, as all disjunct sets and all supersets are counted twice. For this special case, the number of compatible columns is simply $2^n$.


How to generate a consensus sequence from a multi-reference BAM file? I am trying to generate a consensus sequence from a BAM file that was generated by mapping reads to a reference FASTA containing multiple sequences.
Usually, I generate consensus sequences from BAM files using samtools and bcftools:
samtools mpileup -vf reference.fasta filename.sorted.bam | bcftools call -m -O z - > filename.vcf.gz

bcftools index filename.vcf.gz

bcftools consensus -f reference.fasta filename.vcf.gz > filename.consensus.fasta

However, this doesn't work when the reference FASTA contains more than 1 sequence.
Is there a way to generate a consensus sequence from a multi-reference BAM file?
It would also be awesome if the method filled in gaps in coverage with Ns or Xs in the generated consensus sequence but this isn't essential.
Thank you for your help!

This post on biostars has several suggestions.
THe one that seems most straightforward, starting where you are, is to use vcftools:
cat ref.fa | vcf-consensus file.vcf.gz > out.fa

Based on that I think that you should also be able to solve this problem using alfred and GATK.
If all of these suffer the same problem of needing a single sequence (which I very much doubt that they do), you could split the bam / ref files into multiple files, run whatever consensus workflow, and then combine them again. But that doesn't seem like it should be necessary.
See also this post and this post.


Get protein names corresponding to PDB ID I have a list of about 4000 PDB IDs and would like to get the actual names of the proteins (e.g. lactate dehydrogenase, cytochrome c). I tried the batch header section at the Protein Databank Download page but it refused to accept my PDB IDs in formats (xxxx or xxxx.pdb, individually or space-separated) that worked in an interactive search for the protein structure.
Any suggestions?

You can use one of the UniProt Protein APIs.
As you said you have your pdb entries in a text file line by line you can, like this example.txt containing:
1brr
4lzm
2dyi

Using the commandline, you can use a little script like this to download the name, if it is available for the given pdb entry.
while read line;
do
 curl -X GET --header 'Accept:application/json' "https://www.ebi.ac.uk/proteins/api/proteins/pdb:$line" |
 jq -r '.[].protein.recommendedName.fullName.value'  |
 sed "s/^/$line\t/" >> pdb_names.txt;
done < example.txt;

You need to have curl, sed and jq installed on your system.
This gives you following output in pdb_names.txt
1brr    Bacteriorhodopsin
4lzm    Endolysin
2dyi    Ribosome maturation factor RimM

Update
if you want to speed it up, you can run it with parallel.
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" |   jq -r ".[]. .protein.recommendedName.fullName.value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

With the -j option you call how many jobs should run in parallel. The limit of the UniProt API is 200 request per second per user.

Update 7. Nov 2020
To get another info beside the protein name, you need to know how the JSON response from UniProt looks like.
To get also the scientific name, you can run following command:
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" | jq -r ".[] | .protein.recommendedName.fullName.value + \" - \" + .organism.names[0].value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

As result you get this:
1brr    Bacteriorhodopsin - Halobacterium salinarum (strain ATCC 700922 / JCM 11081 / NRC-1)
4lzm    Endolysin - Enterobacteria phage T4
2dyi    Ribosome maturation factor RimM - Thermus thermophilus (strain HB8 / ATCC 27634 / DSM 579)


How to map reads shorter than 32bp with minimap2? I mapped fasta sequences to a reference. The fasta sequences range in length from 17bp up to several hundred bp. I used minimap2 with the following command: minimap2 -ax sr ref.fa seqs.fa, et cetera. Only the sequences greater than 32bp len mapped though.

samtools view -F 4 filename.sorted.bam | \
  cut -f10 | awk '{print length($1)}' | \
  sort | uniq -c | sort -k2 -n

   2304 32
  34931 33
   5443 34
   5343 35
   8811 36
   5485 37
   5417 38

Is it possible to make minimap2 map sequences shorter than 32bp? Do I need to use bwa aln or some other mapper?

Minimap2 doesn't work well with 32bp reads by design. Older mappers like bowtie1 and bwa-aln will be better.


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

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

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

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

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

No, downsampling in general has no use as it reduces power. Peak callers such as macs2 will internally scale data to account for sequencing depth differences.
The IgG itself is only useful for peak calling as it somewhat corrects for the unspecificity of the antibody to precipitate any kind of protein rather than the protein of interest. I hear people say, and personally agree with this by experience with my own experiments, that of course IgG pulls down notably fewer DNA that the target antibody, therefore the library requires notably PCR to produce enough material for sequencing. And this then makes it questionable how representative IgG is as a control. I personally prefer chromatin input to account for the intrinsic sequence/GC bias of the underlying genomic region in the sequencing/mapping. To my knowledge there are currently few established tools to really take these controls into account for downstream analysis such as differential testing via interaction terms because the library composition is simply so different from the IPs making normalization difficult.
Your focus imho should be rather to properly normalize your samples to correct for differences is signal/noise ratio, see for example my post over at biostars on how one could do that beyond naive per-million scaling.


Calculating Isoelectric point from a multifasta file I need to calculate Isoelectric point from a multifasta file, is there any python code or web tool that allows me to do that?

You can use EMBOSS Pepstats for this. It takes a multi-fasta as input and produces a table that contains various statistics related to the protein, including isoelectric point.


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


Error in formation of correlation matrix-(corrupt matrix -- dims not not match length) I am working on a script to create a correlation matrix and mutual_rank matrix from RNAseq gene counts for 57992 genes in 7027 samples. I have already tested this script for a smaller file (for 10,000 genes) and it is working, however, when I am running this script on my final input file, it is giving the following error:

 corrupt matrix -- dims not not match length (this error comes after gene to
 gene correlation step)

Script which i am running is:

E <- read.table("mdat_human_test.txt", header = TRUE, row.names=1)

# gene-to-gene correlation
C <- cor (t(E))
write.table(C, "human_correlation_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

# correlation rank matrix
R <- t(apply(1-C, 1, rank))
write.table(R, "human_rank_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

# mutual rank matrix
M <- sqrt(R * t(R))
write.table(M, "human_mutual_rank_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

q()

You can't compute a correlation matrix when you have more than ~46000 rows, since a standard R matrix can have a maximum of 2^31-1 values. Have a look at packages like bigcor. Alternatively, consider if having 3+ billion correlation values is really useful.


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?

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.


Is the optional SAM NM field strictly computable from the MD and CIGAR? From SAM Optional Fields Specification the NM field is

  Edit distance to the reference, including ambiguous bases but excluding clipping

Assuming both the MD and CIGAR are present, is the edit distance simply the number of characters [A-Z] appearing in the MD field plus the number of bases inserted (xI, if any) from the CIGAR string? Are there any other complications?

Assuming both the MD and CIGAR are present and correct, then yes, you can parse both to get the edit distance (NM auxiliary tag). One big caveat to this is that there's a reason that the samtools calmd command exists, since it's historically been the case that not all aligners have output correct MD strings. It's rare for the CIGAR string to be wrong and that'd be more of a catastrophic error on the part of an aligner. For what it's worth, if the NM auxiliary is absent on a given alignment but present on others produced by the same aligner then it's fair to assume NM:i:0 for a given alignment by default (many aligners only produce NM:i:XXX if the edit distance is at least 1).


Can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence? I am trying to perform de novo structure prediction of a 1135 long protein sequence; can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence?

For DMPfold you might like to look at Figure 5B of the paper, which implies that performance starts to drop off after ~600 residues. Only proteins up to 500 residues were used in the training set, for example.

A protein of that length is almost certainly multi-domain. I would consider doing domain prediction; we have a tool DomPred, but many others are available. Then run the sequences of the domains individually and recombine them later.

I don't know about I-Tasser's ability to do this, perhaps see if they mention it in their papers.