How to input the opposite strand of DNA? I am trying to write a program that will count the repeats of a given DNA strand:
I am using a Counter library, something like this code shown below
from collections import Counter

        s="XXXXXAGCCTGCCAAGCAAACTTCACTGGAGXXXXXTGTGCGTAGCATGCTAGTAACTGCAXXXXXTCTGAATCTTTCAGCTGCTTGTXXXXXTGGGCCTCTCACAAGGCAGAGTGTCTTCATGGGACTTTGATATTTATTTTTGTACAACCTAAGAGGAACAAATCCTTTGACACTGACAAATTGGCTTCCATATTTTATACCTTAATCATCTCCATGTTGAATTCATTGATCAACAGTTTAAGAAAAAAAGATGTAAAAATGCTTTTAGAAAGAGAGGCAAAGTTATGCACAATAACTTCTCATGAAGTCACAGTTTGTTAAAAGTTGCCTTAGTTCACAATAAATAATTATGTATGCTCTATAATTTCAGTGA"

def chunk_string(s, n):
    return [s[i:i+n] for i in range(len(s)-n+1)]

counter = Counter(chunk_string(s, 5))

Repeat supposed to be "declared" as an identical sequence on the same strand (like XXXXX here that occurs 4 times). However, I can also have this repeat (for ex. XXXXX in this particular test case) on the OPPOSITE strand.
How am I suppose to intake the opposite strand in case of using .fasta files? As I know for today - fasta is presented only in a single strand way. Or am I completely mixed the terms?

Yes, you're going to want to take the reverse complement of your string in order to analyze the other strand, that is the standard.
Also, it is worth pointing out that you may want to use a generator expression instead of a list expression for your chunk_string function:
def chunk_string(s, n):
    return (s[i:i+n] for i in range(len(s)-n+1))

This will result in a generator that yields kmers as the Counter constructor iterates through them, which means that you don't need to first construct and store the entire list of kmers first and then pass to Counter but instead can feed them directly to it. Not essential for small strings, of course, but a good practice and useful for larger strings. See this SO post for a better explanation.
You can also use Biopython's reverse_complement function in order to easily obtain the reverse complement.


Why are my Chi-squared test results different from those in a published table? I recently read the paper “A novel long non-coding RNA linc-ZNF469-3 promotes lung metastasis through miR-574-5p-ZEB1 axis in triple negative breast cancer”. In this I see Table1 showing correlation of Linc-ZNF469-3 with different features in TNBC and non-TNBC patients. p-values are calculated based on Pearson's Chi-sq test.

For non-TNBC patients: I calculated the p-value with the Chi-sq test just to have a look if everything is right.

Recurrence:

data <- data.frame(x= c(40,62), y=c(56,65))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 1.1269, df = 1, p-value = 0.2884

Brain:

data <- data.frame(x= c(6,99), y=c(6,119))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 0.096461, df = 1, p-value = 0.7561

You can see that the Recurrence p-value doesn't match with what I see in the table. Others, like Brain, Bone, Liver, Lung match with the p-values I calculated.

And for TNBC patients, not a single p-value matches with my calculations. I don't even get Lung and Bone as having significant p-values with my calculations.

I did it for Lung as shown here:

data <- data.frame(x= c(25,143), y=c(9,35))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 0.8044, df = 1, p-value = 0.3698

But in the paper they mentioned:

  We found that high linc-ZNF469-3 expression was significantly correlated to lung metastasis in TNBC patients (Table 1).

Can anyone tell me what could be the wrong here? Is it my calculation wrong? or in the paper?

Your calculations seem right, perhaps there was an error on their side.

I also looked at the number of samples reported, but they use the same amount of samples in each case. Because they are testing whether Linc-ZNF469-3 is different in several tissues, they should have used a multiple test correction for the p-value (although reporting the original values is usually also ok).

To answer your question in the comments. You can't correlate a categorical variable with a continuous  variable. Or they did another test or they are not explaining clearly what they did.


Is nanopolish worth it since faster polishing software is available? For Oxford Nanopore contigs produced by any long-read assembler has anyone performed any benchmarks to compare the polishing tools racon, nanopolish, or something like wtpoa-cns?

The compute time of nanopolish is quite long and I would rather use another tool if it will suffice. After long-read polishing I will polish more with two rounds of pilon using 200x coverage Illumina shotgun data from the same individual.

Nanopolish is necessary if you want to get high-quality consensus. racon etc don't use signal data. They can't achieve high quality. pilon at times doesn't work well. Even if it worked perfectly, it wouldn't help repetitive regions inaccessible by short reads.


How to verify if the stationary states of a biological model are correct? I'm mathematically modeling yellow fever, and yet, after assembling the model, I need to find the stationary states of the model. I had to do it by hand because wolfram software could not get the necessary values.

My question is, how do I check if my stationary states are in line with my biological model?

I was told that there is a way to verify this by performing a substitution of the s values of the points within my own template code. But I have no idea how to do it.

The code used is taking forever to get back to me as this data behaves in my model

function ypsir = ypsirtry(t,y)
betah = .01;
muh = 0.01;
sigmah = 0.05;
alphah = 0.06;
deltah = 0.05;
mui = 0.07;
mum = 0.07;
betai = 0.03;
lambdai = 0.06;
lambdah = 0.5;
k = 500;

ypsir = zeros(5,1);
ypsir(1) = sigmah*y(3) - muh*y(1) - betah*y(1)*y(5) + lambdah*(y(2)+ y(3) +
y(1));
ypsir(2) = -deltah*y(2) - y(2)*(alphah + muh) + betah*y(5)*y(1);
ypsir(3) = deltah*y(2) - sigmah*y(3);
ypsir(4) = -mui*y(4) - mum*y(4) - betai*y(2)*y(4) + lambdai*(y(4) + y(5))*(1
- (y(4) + y(5))/k);
ypsir(5) = betai*y(2)*y(4) - mui*y(5) - mum*y(5);
end

Along with this calling program

clear;
to = 0;
tf =50;
yo = [99 1 0 1 3];
[t y] = ode45(@(t,y)ypsirtry(t,y),[to tf],yo);
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5))
title('First Model of Yellow Fever')
xlabel('time')
ylabel('susceptible, infected, recovered,Mosquito Non-carrier, Carrier
mosquito')

My stationary states are:

S = (1/betai)*[(muh+alphah)+(sigmah/deltah)*deltah]*(mui+mum)/betai*N*I);

I = sigmah*R/deltah;

R = ((muh/betah + lambdah)*(muh+alphah))/(sigmah+lamdah)-(muh+alphah)*
(deltah/sigmah*deltah)-(betah/muh - lambdah)*(betai*N*I)*
(delta/sigmah*deltah));

N = -k/lambdai((-lamdai+(lambdai/k))*P/-betai*I + lambdai - (2k/lambdai*P)-
mui-mum);

P = betai*N*I/(mui+mum);

At where
S = susceptible
I = infected
R = retweeted
N = non-carrier mosquito
P = carrier mosquito

After much analysis of the model, I came to the conclusion that I will not be able to verify how the equilibrium points behave graphically. Since you will have to create a phase diagram in R ^ 5, there are five equations, and I can not do two-to-two comparison in these cases, because I will always have a function that depends on a variable that will not be placed, but will be in the equation plotted anyway.

But in any case, if it were in three dimensions, for the purpose of information for those who are interested, the verification of the points would be to place the points (which in calculating have both numbers and not functions in function of other functions), within the model created, hence plot the phase diagram and check whether the arrows converge or not to the point. If they are converging, they are correct.

Regarding the model is taking forever, I took the logistics part and ran well!


Output when compute GC content is different between using user input and using string I am learning bioinformatics and i was solving a problem on rosalind.info.
Input is a string consist of multiple DNA strand names and the DNA sequences, and Output is the name of the sequence that has the highest GC percentage and the GC percentage.
this is the input:

>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

this is the correct output:

Rosalind_0808
60.919540

This is my code (sorry for this filthy mess):

DNA = input()
wholeRosalind_DNA = DNA.splitlines()
Rosalind = []
Rosalind_DNA = []
for i in wholeRosalind_DNA:
    if "Rosalind" in i:
        Rosalind.append(i)
Countstring = []
maxGCcontent = []
GCcontent = []
def eachRosalind_DNA(DNA):
    for i in range(0,len(Rosalind)):
        if Rosalind[i] in wholeRosalind_DNA:
            a = wholeRosalind_DNA.index(Rosalind[i])
            if a == wholeRosalind_DNA.index(Rosalind[-1]):
                b = len(wholeRosalind_DNA)+1
            else:
                b = wholeRosalind_DNA.index(Rosalind[i+1])
            Rosalind_DNA.append(wholeRosalind_DNA[a:b])
    for i in Rosalind_DNA:
        Countstring.append(''.join(i))
    for i in Countstring:
        bases = len(i)-14
        CG = i.count("C") + i.count("G")
        GCcontent.append([i[1:14],(100*(CG/bases))])
    for i in GCcontent:
        maxGCcontent.append(i[1])
    result = max(maxGCcontent)
    for i in GCcontent:
        if i[1] == result:
            return (i[0]),(round(i[1],6))
print(eachRosalind_DNA(DNA)[0])
print(eachRosalind_DNA(DNA)[1])

but i get this instead:

  File "GC_Content.py", line 26, in eachRosalind_DNA
    GCcontent.append([i[1:14],(100*(CG/bases))])
ZeroDivisionError: division by zero

i will get the desired output if i just paste the input into the variabe "DNA" without using user input. Please help me fix the code and distinguish those two ways of getting input

Is it cheating to use Biopython(?), will a solution might be helpful anyway. This assumes you have your sequences stored in the seq.fa FASTA file:

from Bio import SeqIO
from Bio.SeqUtils import GC

print(max(((r.id, GC(r.seq)) for r in SeqIO.parse("seq.fa", "fasta")),
          key=lambda x: x[1]))

Output:

('Rosalind_0808', 60.91954022988506)


What is this 5-column sequencing format from 2009? There is high throughput sequencing data here, and I don't know what format it is in.

It was submitted in 2009, and the description says the following:

Library strategy:     ncRNA-Seq
Library source:   transcriptomic
Library selection:    size fractionation
Instrument model:     Illumina Genome Analyzer II
Description:  CIPPNK, tar file of Illumina *_seq.txt files provided as supplementary file

I got the archive here:

ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM455nnn/GSM455387/suppl/GSM455387%5FWT%5FCIPPNK%5Fseq%5Fs1%2Etar%2Egz

Inside, there are 330 files from s_1_0001_seq.txt to s_1_0330_seq.txt which are tab-delimited text files where the first column is always 1, the second has the number found in the file name, then 2 mysterious integers, and then what looks like a read of length 36, with sometimes a dot instead of a letter:

$ head s_1_0330_seq.txt
1   330 690 785 TTCCTACATTGTTCCCCCATGCTGTTGGCACCATCA
1   330 44  145 TTTTTATCACGAGTTTTAAATCTGTAGTCACCATCA
1   330 53  141 AATAATGCATAACAAAACGGAATCTGTAGAA.AAA.
1   330 784 461 TAATTGTAGTGATTGATCAATCTGTAGGCACCATCA
1   330 588 634 TATTATGCACATTTTCTAGTTCACTGTAGGCACCAT
1   330 718 678 TTACATGTTTCGGGTAGGAGCCTGTAGGCACCATCA
1   330 635 834 TGTGATCATTAGTTCAAAGCCCCCTGTCGGCACCCT
1   330 494 523 TGAAAATCAAAAATGCTGAACTGTAGGCACCATCAA
1   330 393 783 TTTTTTTTTAAATTTAAAAAAACTGTAGGCACCATC
1   330 48  148 GTTTAACCGTGTAGACGTTGGTTTCTGTAGGCACCA

This seems to be some "old" high-throughput sequencing format. I think someone in a 2008 message in sequanswers was dealing with this type of file:

http://seqanswers.com/forums/showpost.php?p=1841&postcount=8

What is this format that seemed so standard back then, that the authors did not give more info than describing the files as "Illumina *_seq.txt files"? I do not dare ask them for such a trivial question (the given contact is a Nobel prize laureate and is probably too busy to answer random bioinformatics questions).

In particular, what are the columns 3 and 4, and what do the dots mean?

This is an early Solexa/Illumina sequencer format. The columns are the identifier to location on the flowcell.  I believe the "." was the original placeholder for an unread base, which has been replaced by and "N" in current Illumina sequencing output.
From on http://www.crg.eu/en/content/processing-and-analysis-illumina-sequencing-data.

seq.txt (Gerald)
4 1 23 1566 ACCGCTCTCGTGCTCGTCGCTGCGTTGAGGCTTGCG
The first 4 columns are an ID that provides location details of the cluster in the flow cell, followed by the read sequence. Seq.txt file size typically ranges between 250 Megabytes and 1 Gigabyte. One file is generated per lane.


Does read sequence identifier need to be unique for proper alignment? I am preprocessing a set of scTHSSeq reads in which each read sequence identifier has been replaced with the cell barcode, e.g.:

@CGAGGCTG_TATCCTCT_AGTCCA
...
...
...
@CGAGGCTG_TATCCTCT_GGTACG
...

The read sequence identifier has been replaced with the cell barcode to be able to demultiplex alignments after bwa into a different bamfile for each cell.

I ran the alignment using bwa mem and I received no errors. However I would like to know if there are caveats or reasons why this approach may lead to misleading results.

I will then proceed to remove PCR duplicates with samtools rmdup executed in each cell's bamfile and then I'd run peak calling using spp.

What tool did you use to replace the read names? Can you modify that to split up the reads into separate files or encode the cell barcode in the read description instead of the read name?

Typically is it assumed by many tools that the FASTQ read names are unique. you could have issues with BWA and split alignments but I don't know if that pertains to your data.  I think the best you can do is the sanity check the results at each step of the process.

In the fastq format (https://en.wikipedia.org/wiki/FASTQ_format) on the first line denoted by the '@' there can be two values separated by a space.  The first is the name and the second is the description.

So an example would be:

@Original_read_name CGAGGCTG_TATCCTCT_AGTCCA
ATCGATCATGATGACA
+
@&%*@(%@%*($!##!


How does Li and Durbin's BWA paper compare alignment programs on real data? Li and Durbin's "Fast and accurate short read alignment with burrows-whleeler transform" found here, says:

  We evaluate the performance of BWA on ... real paired-end data by
  checking the fraction of reads mapped in consistent pairs and by
  counting misaligned reads mapped against a hybrid genome.

The data for the real genome aligned to is in table 2.

Which pairs are referred to here? Are these the different software tools aligning the same corresponding reads meaning that they aligned the same sequences to similar locations? Are these paired ends? Or is this something different?

What they are pairing is the RNAseq data with the reference genome: as the legend of the table 2 says:

  The 12.2 million read pairs were mapped to the human genome

These softwares align the same read to the same genome. They aren't comparing if it is mapped to the same location, only how many mismatches has each program.


Import a tab-separated file with differing numbers of elements in each row; prokka output I am using prokka to annotate a bacterial genome:

prokka ecoli.fa

Prokka is outputting a tab-separated file (called PROKKA_12142017.tsv) with differing numbers of elements in each row:

locus_tag    ftype    gene    EC_number    product
BOFHCHGE_00001    CDS    hypothetical protein
BOFHCHGE_00002    CDS    rdgB_1    3.6.1.9    dITP/XTP pyrophosphatase
BOFHCHGE_00003    CDS    hypothetical protein
BOFHCHGE_00004    CDS    hypothetical protein
BOFHCHGE_00005    CDS    hypothetical protein

Does anybody know how I can get round this to import into a dataframe in R?

Is there a way to ask prokka to add something in the gene column?

You don't really need to get around anything, your example can be loaded correctly with read.delim(). The rows with missing values are filled in with blanks. You may, however, prefer the readr package, which handles this a bit more elegantly (it'll tell you which values were missing and not do the annoying string->factor conversion).


Calculate molecular mass from chemical formula C9H15N4O8P with Python? Can I calculate the molecular mass from a formula like this C9H15N4O8P with Python e.g. using rdkit? I think rdkit.Chem.Descriptors.MolWt needs a molecule object. Though, to generate that I would need e.g. a SMILES string. Is there any other python module that can do the trick?

I found this: http://verahill.blogspot.com/2012/10/molecular-weight-calculator-in-python.html

As you mentioned, rdKit does not accept a chemical formula as input. The main reason for this are probably isomers. You could however use rdKit to initiate all the atoms, get their mass and add them together.

import re
from rdkit import Chem

def get_mass(formula):

    parts = re.findall("[A-Z][a-z]?|[0-9]+", formula)
    mass = 0

    for index in range(len(parts)):
        if parts[index].isnumeric():
            continue

        atom = Chem.Atom(parts[index])
        multiplier = int(parts[index + 1]) if len(parts) > index + 1 and parts[index + 1].isnumeric() else 1
        mass += atom.GetMass() * multiplier

    return mass

Then...

>>> get_mass("C9H15N4O8P")
338.21299999999997


Differences between NetAffx Hg-U133 Plus 2.0 Annotation file versions From what I was able to observe, the annotation files (CSV) for NetAffx Hg-U133 Plus 2.0 is updated regularly.

However, from glancing at the CSV files (NA35: 2015 & NA22: 2007), I can't seem to notice any glaring difference between the two.

Does anyone know what information is actually being updated? I assume it may be gene identifiers for different databases, but wanted to make sure.

Thanks.

Below are the files I downloaded and examined from Affymetrix website (you need an account to download archived annotation files)

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

After checking the files myself, it does seem like gene identifiers indeed have been updated. Some that were missing have been added with identifiers and some have additional identifiers now.

Among some other things that were updated from 2007 include:

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

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


Are there databases or resources that contain expression of proteins or genes before and after cancer treatment? Are there databases that contain information about gene or protein expression before and after targeted therapy or chemotherapy?

I am curious about different cancer cell lines or cancer types.

Any suggestion on a set of resources I should look into?

One such database/tool is “[PharmacoGx][1]: an R package for analysis of large pharmacogenomic datasets.” Bioinformatics (Oxford, England).

They brought together a large compendium of experimental data (cell lines and the like if I remember correctly from a talk), standardized these and made available within a special object class. I think there is genomic and transcriptomic data, not sure if there is proteomic data though.


How are PDO and PDX used in computational and predicative models for tumour biology? There are wet methods: Patient-derived models: Patient Derived Xenograft (PDX) and Patient Derived Organoids (PDO) to reflect tumor biology.

Are there any databases or computational tools that use the outcomes from PDO/PDX experiments to create a predictive computational model for cancer or other diseases?

Potential pitfall!

I'm not sure about predictive models, but you need to be aware of a potential pitfall in blindly aligning PDX or PDO based sequencing data without first removing contaminating host organism reads, as otherwise these will lead to a lot of false positive variants caused by miss alignment.  In my experience even a small mount of host material can lead to a ten fold increase in called variants due to miss aligned host reads looking like true variants.  I'd recommend using Xenome, source.

Note Xenome was designed for DNA sequencing it's not splice-aware, so for RNA-Seq there is no optimal solution other than removing reads which align to the host organisms genome.  Although the issue here is that in conjunction with a splice-aware read aligner synteny between the host and grafted genomes might create some interesting problems.  However for RNA-Seq provided FACS or similar confirms low-levels of host contamination I expect levels of expression will not be badly affected.  Although this really needs investigating.

Finally and rather annoyingly Xenome produces none-standard FASTQ so you'll need to fix it's output with:

awk '{if (NR % 4 == 1) print \"@\"$0; else if (NR % 4 == 3) print \"+\"$0; else print $0 }' as reported here on seqanswers.


Merging / Concatentating BLAST Databases I was wondering if any one has had success merging BLAST databases from multiple references?

For a concrete example I have hg19 and a bunch of viral and bacterial genomes that all have been index separately but I want to do a single alignment to all of them.  Is there a way I could say combine hg19, virus_A, virus_D and bacteria_X into a single database without having to regenerate?  Since the problem could get combinatorial quickly and making a db with hg19 is time consuming I was wondering if there was a way around it.

This is similar to another ticket but asking for BLAST instead of alignment indexes.

There's an app for that, apparently.

  Often we need to search multiple databases together or wish to search a specific subset of sequences within an existing database. At the BLAST search level, we can provide multiple database names to the “-db” parameter, or to provide a GI file specifying the desired subset to the “-gilist” parameter. However for these types of searches, a more convenient way to conduct them is by creating virtual BLAST databases for these. Note: When combining BLAST databases, all the databases must be of the same molecule type. The following examples assume that the two databases as well as the GI file are in the current working directory.

  Aggregate existing BLAST databases

  To combine the two nematode nucleotide databases, named “nematode_mrna” and “nematode_genomic", we use the following command line:

$ blastdb_aliastool -dblist "nematode_mrna nematode_genomic" -dbtype nucl \
  -out nematode_all -title "Nematode RefSeq mRNA + Genomic

So, you can use the blastdb_aliastool to create merged databases just like you asked for.


Filtering Sequences (entries) by headers ID from a Fasta file database First of all, sorry if this question has been posted previously (I could not found a solution accessing the previous Q&A).
I have a fasta file as follow:
>LNIV02000036.519060.520603 Pseudomonas aeruginosa
GAACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAA...

>LDLH01000045.361.1876 Pseudomonas aeruginosa
AGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATG...

>LNMQ01000087.5825.7394 Enterococcus faecium
TTTTTATGAGAGTTTGATCCTGGCTCAGGACGAACGCTGGCGGCGTGC...

...

I would like to filter sequences or extract entries for a specific species, let's say Pseudomonas (I have others in mind...). I tried a few things but I think the main problem is my skills associated with bioinformatics. I would like to create a file content Pseudomonas entries (~9k entries in the current database). Any suggestion or tip?
Thank you for your time.
Best,
Fabricio

Using Biopython:
from Bio import SeqIO

for record in SeqIO.parse("input.fa", "fasta"):
    if record.description.split()[-2:] == ["Pseudomonas", "aeruginosa"]:
        print(record.format("fasta"))


Generate table for total number of SV events per sample I have called and annotated structural variants and I have a table by name of samples and structural events like below

Sample  Type
t_005   DEL
t_005   INS
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DUP
t_005   DUP
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DUP
t_005   BND
t_005   DEL
t_005   DEL
t_005   BND
t_005   DEL
t_005   DEL
t_005   DUP
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL

How can I count the number of events to get the below table?

> head(SV_Counts[,1:4])
               t_005      o_005       t_013         o_013
    DEL          27          40          25          30
    DUP          11          10          13          22
    BND          23          20          10          12
    TRA           3           2           4           4
    >

I tried this from google but I failed to adapt my data to that

    > melt(table(a))
   Sample Type value
1   o_036  BND   280
2   o_005  BND     8
3   o_013  BND   160
4   o_021  BND   144
5   o_036  BND     0
6   t_005  BND    15
7   t_013  BND   108
8   o_036  DEL    63
9   o_005  DEL    15
10  o_013  DEL   121
11  o_021  DEL   243
12  o_036  DEL     0
13  t_005  DEL    40
14  t_013  DEL    76
15  o_036  DUP    34
16  o_005  DUP     4
17  o_013  DUP    57
18  o_021  DUP    88
19  o_036  DUP     0
20  t_005  DUP    14
21  t_013  DUP    48
22  o_036  INS    58
23  o_005  INS     0
24  o_013  INS     4
25  o_021  INS    10
26  o_036  INS     1
27  t_005  INS     1
28  t_013  INS     2

But hereafter I don't know what to do

So you have data in "long format" and want to convert it to "wide format". You can use spread() from tidyr for that:

require(tidyr)
# Assuming your input is called `foo`
spread(foo, Sample, value)


How to know if ID correspond to nucleotide sequence or protein with a script I have a text file that contains a list of IDs (314 sequences):
AVP78031.1
AVP78042.1
ATO98108.1
ATO98120.1
ATO98132.1
...

My goal is to make a script (maybe using Python or Perl) to check in the list if all the IDs are nucleotide or protein sequences.
For example:
AVP78031.1 -> protein (this is a nucleotide sequence, I change nucleotide for protein to show an example).
AVP78042.1 -> nucleotide
ATO98108.1 -> nucleotide
ATO98120.1 -> nucleotide
ATO98132.1 -> nucleotide

Any idea to do a script?
Thank everybody!

If these are all GenBank or RefSeq accessions, you can use Entrez Direct for this as shown below:
$ cat accs.txt
ATO98108.1
ATO98120.1
ATO98132.1
AVP78031.1
AVP78042.1
$ cat accs.txt | epost -db nuccore | efetch -format acc
## no output because none of them are nucleotide accessions
$ cat accs.txt | epost -db protein -format acc | efetch -format acc
AVP78042.1
AVP78031.1
ATO98132.1
ATO98120.1
ATO98108.1

NOTE: This will work only if the accessions are currently live because epost does not find any suppressed accessions. For example:
$ cat accs.txt
NM_002826.3
NM_002826.4
NM_002826.5
$ cat accs.txt | epost -db nuccore -format acc | efetch -format acc
NM_002826.5

Here, all three accessions are valid nucleotide accessions but only the last one, NM_002826.5, is alive.
An alternate way is to use the accession prefixes defined here and come up with an appropriate regular expression query.


SC-RNA seq percent.mito Is that possible that the percent.mito level is zero? I have done vlnplot and it shows that its percent.mito level is zero. Does that mean my  mitochondrial gene didn't detect?

It is violin plot. And the command is:

I've just come across a similar issue. It seems the mouse and human mitochondrial genes use mixed capitalizations. See the following PR for details:

https://github.com/Hoohm/dropSeqPipe/pull/43

Bottom line is, if you ignore the case in the grep pattern it might work (assuming you use human/mouse):

mito.genes <- grep("^mt-", rownames(shox2crepos_mtmg@data), value=TRUE, ignore.case=TRUE)


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

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

solution from @sören:
I came up with this solution for my problem:
from Bio.KEGG import REST as kegg

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

_get_kegg_v = np.vectorize(_get_kegg)

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


How to obtain .bed file with coordinates of all genes I want to get a .bed file with the genes' names and canonical coordinates, also I would like to have coordinates of exons, too. I can get the list from UCSC, however, if I choose UCSC Genes - knownCanonical, I can not extract coordinates of exons. If I use other options - I am getting coordinates of as many transcriptional isoforms as were detected while I need only one canonical form.

How can I get such BED file?

Via Gencode and BEDOPS convert2bed:

$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

You can modify the awk statement to get exons, by replacing gene with exon.

BEDOPS: https://github.com/bedops/bedops

This is based off an answer I wrote on Biostars, which includes a Perl script for generating a BED file of introns from gene and exon annotations: https://www.biostars.org/p/124515/#124522


How do I add a colour annotation bar to the heatmap generated by DoHeatmap function of Seurat v2? I am using Seurat v2 for professional reasons (I am aware of the availablity of Seurat v3). I am clustering and analysing single cell RNA seq data. How do I add a coloured annotation bar to the heatmap generated by the DoHeatmap function from Seurat v2? I want to be able to demarcate my cluster numbers on the heatmap over a coloured annotation bar. I want to have the an image like this:

 (taken from Satija V3 tutorial).

My original code:

library(Seurat)
DoHeatmap(object = object, genes.use = top20_pos$gene)

produces the following image which is not what I want:

I am aware that DoHeatmap function returns a ggplot2 object but after searching the internet for a long time, I have not been able to find  which geom layer I can add to achive my objective.

On the other hand, I have tried to generate a separate coloured annotation bar from another package called ComplexHeatmap with the intention of cropping it to fit the cluster demarcations on my heatmap generated with DoHeatmap but this has proven to be very time-consuming and inaccurrate:

library(ComplexHeatmap)

df = data.frame(type = c(rep("a", 18), rep("b", 3), rep("c", 2), rep("d", 1), rep("e", 1)))
ha = HeatmapAnnotation(df = df, col = list(type = c("a" =  "blue", "b" = "red", "c" = "green", "d" = "cyan", "e" = "magenta")))
draw(ha, 1:25)

which produces:

I will be very grateful for any hint on how to achieve my goal automatically and accurately.

Thank you in advance for your help.

I am not sure if it's a v3 or also v2 thing, but have you tried by setting group.by and label?

From ?Seurat::DoHeatmap:

group.by   A vector of variables to group cells by; pass 'ident' to group by cell identity classes
label      Label the cell identies above the color bar

If you prefer external tools, you can use pheatmap and setting annotation_col e annotation_colors.

First define a data.frame (df) containing the data you want in the bar (clusters in your case)

Description Condition
liver       Tumour
lung        Normal

Then you have to assign the colours

anno <- as.data.frame(df)
anno_colori <- list(
   Condition = c(Normal = '#FFC000', Tumour = '#0070C0'),
   Description = c(liver = '#D55E00', lung = '#66CC99')
)

pheatmap(mat,
     annotation_col = anno,
     annotation_colors = anno_colori,
)

The code above will add 2 bars on the top of your heatmap, one for Condition and another for Description.


Is there a "definitive" database for Mobile Genetic Elements? This isn't quite my speciality, but I've been tasked with finding any evidence of Mobile Genetic Elements (like transposable elements, retrotransposons, bacteriophages, etc.) within some WGS samples.

Is there a "definitive" database used with sequences? How "species-specific" would this be?

Furthermore, are there standard algorithms used to find these MGEs? I suspect we don't have a comprehensive database of all possible MGEs, but it may be possible to query against the sequences we do indeed have (or detect any repetitive elements within the genome).

Would advise doing Denovo identification method structure-based. Then use these public databases (  RepBase, Dfam ) to look for high similarity or homology exploration, this possibility if your species, the repeats not yet discovered, or deposited in the databases.

If the species well-investigated earlier and the mobile element recognised, I assume it will be a straightforward task using species-specific search via Repeatmasker tool with the -species flag to obtain the MobileElement annotation in the genome.

Let me know if you need any data from Repbase related to your species, i have an entree.


What happens if a major bug is discovered in a bioinformatic package that has been used in published literature? Yesterday I was debugging some things in R trying to get a popular Flow Cytometry tool to work on our data. After a few hours of digging into the package I discovered that our data was hitting an edge case, and it seems like the algorithm wouldn't work correctly under certain circumstances.

This bug is not neccessarily the most significant issue, but I'm confident that it has occurred for other users before me, perhaps in a more insidious and less obvious way.

Given that since this tool was released, literature has been published that utiltised it, how does the profession handle the discovery of bugs in these widely used packages?

I'm sure that given the piecemeal nature of a lot of the libraries out there, this is going to happen in a much more significant way at some point (i.e. invalidating a large number of published results)

For some context I'm a programmer who's been working across various other parts of software development, and is quite new to Bioinformatics.

I prefer to treat software tools and computers in a similar fashion to laboratory equipment, and in some sense biology in general. Biologists are used to unexpected things happening in their experiments, and it's not uncommon for a new discovery to change the way that people look at something. Things break down, cells die off quicker on a Wednesday afternoon, results are inconsistent, and that third reviewer keeps on about doing that thing that's been done a hundred times before without anything surprising happening (just not this time). It's a good idea to record as much as can be thought of that might influence an experiment, and for software that includes any input data or command line options, and especially software version numbers.

In this sense, a discovered software bug can be treated as a new discovery of how the world works. If the discovery is made public, and other people consider that it's important enough, then some people might revisit old research to see if it changes things.

Of course, the nice thing about software is that bugs can be reported back to the creators of programs, and possibly fixed, resulting in an improved version of the software at a later date. If the bug itself doesn't spark interest and the program gets fixed anyway, people unknowingly use newer versions, and there might be a bit more confusion and discussion about why results don't match similar studies carried out before the software change.

If you want a bit of an idea of the biological equivalent of a major software bug, have a look at the barcode index switching issue, or the cell line contamination issue.


how to write the script to know each sample's cluster number? I have followed the Seurat Tutorial:Integrating stimulated vs. control PBMC datasets to learn cell-type specific responses by using my data. And I am confused about how to write the script to calculate the cell number of each sample on each cluster.

I have run TSENPlot to know the samples and I'm using another TSENPlot to get the clusters. If I want to know find the cluster1 of each sample's number how can I write the script?

Assuming you are following the basic seurat guide:

pbmc@ident

You can access the 'ident' slot of your Seurat object and you will get a vector of the cluster IDs and the names will be the respective cells.

head(seuratobj@ident)
Cell1A7 Cell1C4 Cell1D1 Cell1D2 Cell1D3 Cell1D5
0 3 3 1 3 1
Levels: 0 1 2 3 4


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

I totally got 20 modules named with different colors stored in dynamicColors. Like below I first calculated eigenvalues and then used it and trait information to check the association between Modules and trait data.
MEs0 = moduleEigengenes(datExpr, dynamicColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

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

Is Pearson correlation a right way to check the association?

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

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

Answer from @m, converted from comment:
Personally I'd build a tree. Whats your question patients or the relationship between RNAseq 'signatures'? Its not the definitive analysis, but its good mining.
Pearson is doable but how many Pearson's? You'll end up with a 2x2 matrix. How will you resolve the matrix clustering? A maximum likelihood phylogeny solution would work here making use of the binary format. There are other approaches.
Spearman correlation would be better than Pearson's you end up with a 2x2 matrix.
Answer from @peter-langfelder, converted from comment:
You can use Pearson correlation, it is equivalent to carrying out an equal-variance t-test for the module eigengene between the two groups specified by the binary variable.


plotting two heatmaps with the same order of genes I expected the same pattern but here I am not able to compare the patterns as the order of genes does not seem the same. I mean how I can have two heat maps on which the order of genes are the same so I woulds be able to compare the block of yellow and blue colors in two heat maps.
I have used this code
require("RColorBrewer")
myCol <- colorRampPalette(c("dodgerblue", "black", "yellow"))(100)
myBreaks <- seq(-2, 2, length.out=101)
heat <- t(scale(t(sc_DEGG)))
hr <- hclust(as.dist(1-cor(t(y), method="pearson")), method="complete")
heatmap.2(
    heat,
    Rowv=as.dendrogram(hr),
    Colv=as.dendrogram(hc),
    col=myCol,
    breaks=myBreaks,
    main="Title",
    key=T,
    keysize=1.0,
    scale="none",
    density.info="none",
    reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
    trace="none",
    cexRow=0.2,
    cexCol=0.8,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

I  have made a pair of heat maps based on your code in complexheatmap package. The right heat map is lighter than the left one. I don't know which one has been clustered based on the another one. How I can make the right one darker so that both  seem more similar??

ComplexHeatmap is built for plotting side-by-side heat maps with the same clustering - you use the + notation, similar to ggplot2. To use the same example data as @b.nota:
library(ComplexHeatmap)
# First matrix
set.seed(2)
m <- replicate(8, rnorm(26))
rownames(m) <- letters[1:26]

# Second matrix
set.seed(3)
m2 <- replicate(8, rnorm(26))
rownames(m2) <- letters[1:26]

Heatmap(m) + Heatmap(m2)

And the results:

There are lots of examples in the extensive documentation.
EDIT:
Since you have a few genes with high values in the 'T' matrix, they are washing out the colour scheme for the rest of the heatmap. You could set the max to be lower, and accept that those genes will saturate - the resulting heatmap will be less 'honest' about the data, but might capture what you want to see. To change colour scheme requires a call to the colorRamp2 function of the circlize package.
library(circlize)
mycol <- colorRamp2(c(-2,0,2), c("dodgerblue", "black", "yellow"))
Heatmap(m, col=mycol, cluster_columns = FALSE) +
    Heatmap(m2, col=mycol, cluster_columns = FALSE)

This gives:


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

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


How to edit the headers of multiple fasta files from multiple folders? My directories are organized as follows: one main directory, in which I have multiple directories that end with a number ranging from 314 to 727, followed by .3   . For example, 'mgm4761314.3'. Within these directories, I have fasta files named after a similar pattern (for example, mgm4761314.3.299.screen.passed.fna is in the directory mgm4761314.3).

There is only one such file (.fna) per directory.

What I've been trying to do is to write a loop that would go through all these folders and fasta files, and add the short number associated with the folder/file name (eg. 314, 315...727) to the headers. For example, for my first folder/fasta file, instead of:

>seq1
>seq2

I want:

>314_seq1
>314_seq2

The output would be exported to a file called, for instance, 314.fna.

I've sort of made it work for one number at a time, but I have plenty of files, so I'd appreciate if anyone could help me with the loop part! This is what I tried (and it didn't work). Using ubuntu.

for i in {314..727}; do sed "s/>/>${i}_/" [...]/*${i}.3/mgm4761${i}.3.299.screen.passed.fna > [...]/*${i}.3/${i}.fna; done

You were so close! This should do the trick:

for i in {314..727}; do
    sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna
done

Note that I am using -i which makes sed edit the original file. Because I use -i.bak, this will create a backup file with the original data. Once you are sure the command does what you need, you can delete all of these files with rm */*fna.bak.


How to split bedfile into non-overlapping regions and compute aggregation function on duplicate segments I have a bed-like file in which I have several overlapping regions with associated features. I want to split the overlapping regions into disjoint regions in order to process the file (e.g. with pandas), group by regions and compute an aggregate function.

The input data format is the following (overlapping are marked with *):

chr9    8855925    8986425    -1  *
chr9    8855965    8986325     0  *
chr9    8893857    9022105    -1

I want to split into disjoint intervals, then I want to groupby per-interval and compute an aggregation function (e.g. the mean of the last field):

chr9    8855925    8855965    -1
chr9    8855965    8893857  -0.5
chr9    8893857    8986325 -0.67
chr9    8986325    8986425    -1
chr9    8986425    9022105    -1

I was thinking of using bedtools and pandas to aggregate. However, I was unable to find the right tool to do this job in the bedtools/pybedtools documentation.

To make disjoint intervals, you could use BEDOPS bedops --partition, piping to bedmap --mean to get the mean signal over disjoint regions.

Starting with the input bedgraph file, convert it to five-column BED with GNU awk, putting the signal in the fifth column per UCSC convention:

$ awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }' /tmp/in.bedgraph | sort-bed - > /tmp/in.bed
$ cat in.bed
chr9    8855925 8986425 .   -1
chr9    8855965 8986325 .   0
chr9    8893857 9022105 .   -1

The disjoint set from these intervals would look like this:

$ bedops --partition /tmp/in.bed
chr9    8855925 8855965
chr9    8855965 8893857
chr9    8893857 8986325
chr9    8986325 8986425
chr9    8986425 9022105

We can pipe these intervals to bedmap --mean with the original five-column BED file as a "map" file, calculating the mean signal over the disjoint intervals:

$ bedops --partition /tmp/in.bed | bedmap --echo --mean --delim '\t' - /tmp/in.bed > /tmp/answer.bedgraph
$ cat /tmp/answer.bedgraph
chr9    8855925 8855965 -1.000000
chr9    8855965 8893857 -0.500000
chr9    8893857 8986325 -0.666667
chr9    8986325 8986425 -1.000000
chr9    8986425 9022105 -1.000000

The --mean option calculates the arithmetic mean. Other aggregation functions are available to process the signal or score (fifth) column in map BED files, like --stdev, --median, --mad, etc.

When calculating statistics, the default precision is six digits. The --prec <n> option can be used to specify desired precision in the score result.

See bedmap --help or the online documentation for a full description.

To avoid making and cleaning up intermediate files, bash process substitutions could be used to generate a one-liner:

$ bedops --partition <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') | bedmap --echo --mean --delim '\t' - <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') > answer.bedgraph

This can be useful for cleaner pipelines or with processing very large, whole-genome scale inputs on systems with limited disk space.


Is there a tool to convert a variableStep .wig to fixedStep one? I have a .wig file in this format:

variableStep chrom=chr1 span=20
59965   7
59985   10
60005   10
60025   10

and I would like to format it to something like this format:

#bedGraph section chr16:0-103584
chr16   0       59901   0
chr16   59901   59927   1
chr16   59927   59931   2

Please assume, I am working with a reference hg19 genome, in case the info isn't deducible from the variableStep .wig file. I couldn't find anything in WiggleTools

This is a temporary solution I came up with that worked for my specific case:

def convert_line(line, chrom, span, start=None, fixedStep=None, count=None):
        new_line=''
        if fixedStep is None or count is None or start is None:
                line_arr  = line.split()
                start = line_arr[0]
                end   = int(line_arr[0]) + int(span)
                val   = line_arr[1]
                new_line = '{0}\t{1}\t{2}\t{3}\n'.format(chrom, start, end, val)
        else:
                line_arr  = line.split()
                val   = np.float(line_arr[0])
                start = int(start) + int(fixedStep)*int(count)
                end   = start + int(fixedStep)
                new_line = '{0}\t{1}\t{2}\t{3}\n'.format(chrom, start, end, val)

        return new_line

def convert_file(input_file, output_file):
        chrom=None
        span=None
        fixedStep=None
        count=None
        start=None
        with open(input_file, 'r') as fin, open(output_file, 'w') as fout:
                for line in fin:
                        if len(line)<1 or line[0] == '#':
                                continue
                        elif 'variableStep' in line or 'fixedStep' in line:
                                for key in line.split():
                                        if 'chrom=' in key:
                                                chrom = key.split('chrom=')[-1]
                                        if 'span=' in key:
                                                span = key.split('span=')[-1]
                                        if 'step=' in key:
                                                fixedStep = key.split('step=')[-1]
                                        if 'start=' in key:
                                                start = key.split('start=')[-1]
                                                count = 0
                        else:
                                line_out = convert_line(line, chrom, span, start, fixedStep, count)
                                fout.write(line_out)
                                if count is not None:
                                        count += 1

input:

variableStep chrom=chr16 span=20
59965   7
59985   10
60005   10

output:

chr16   59965   59985   7
chr16   59985   60005   10
chr16   60005   60025   10
chr16   60025   60045   10


Problem with classification model of genomic data: every machine learing model predicts wrongly almost always the same subset of dataset First of all, I'd like to apologize for any spelling or grammar mistakes.

I'm having a problem using R for a classification problem. My dataset contains ~300.000 genomic data, and the features are DNA-related features (number of dinucleotides, number of trinucleotides, the CG Content, and some more). In conclusion, I have a dataset of 300.000 rows and 84 columns (columns= features). The 84th feature is basically the classification variable (there are two classes: class 1 and class 2).

I split the dataset in train set and test set and fit the train set to create a model using different methods, like Random Forest, XGBoost, Adaboost, etc. and check the accuracy using the test set.

In every model, the accuracy of the prediction on the test set is always around 86-87%. What I noticed is that, in the test set, for EVERY MODEL, almost the same number of rows (data) is predicted right (the model predict its true class), and always almost the same data of rows (data) is predicted wrong (the model predict wrong class. For example its true class is 1 but the model predicts 2).

If I take the set of data that are predicted right for every model I've created (250.000 data from the 300.000), and build a new model out of them (splitting them again in a new train and test set), I have an accuracy of 98%.

That means that the other 50.000 data (from the 300.000), add an error at every model I'm trying to create.

But I can't just ignore them. I thought that it would be a way to create two differents models, but doesn't know how to figure this out.

Then I thought that there could be a way to find which variables from the two subsets have strong differences, but even if I figure this out, I don't know what's the way for example to scale the original dataset to prevent those differences from keeping the accuracy in 87% in every model (using the original dataset of the 300.000 rows).

In conclusion, I have a dataset that pretends like "two datasets". It consists of two subsets: the first one is always predicted right from any model (machine learning algorithm) I'm trying, and the second one is always predicted wrong. How could I build a strong model taking into account the differences?

I also took a look at the box plots of the features. The most important feature for every model is called "CG" (number of dinucleotide of cytosines and guanines in DNA). In the following image you can see the boxplot of CG feature separating the original dataset in two datasets (T = for the dataset that is predicted right in every model) and (F = for the dataset that is predicted wrong in every model). 1 stands for data of class 1 and 2 stands for data of class 2. So, T1 means data from the right predicted dataset of class 1, T2 means data from the right predicted dataset of classes 2, T means data from the right predicted dataset for both class 1 and 2, etc.

As I'm a new beginner in R and in machine learning, I hope you do understand my problem. Thanks sincerely.

You could be seeing the effect of SNPs on haplotype-driven conditions (i.e. where a combination of multiple variants predicts the class, rather than a single variant).

Do you have any method to determine the age of particular variants (e.g. via an LD plot to find inconsistencies between variants and the surrounding block structure), so that you can see if that correlates with the likelihood of error?


How to quantify similarity of genomes and find differences in set of S aureus genomes? I have around 500 annotated proteomes of different bacterial strains and would like to quantify their similarity (or difference). I found gt genomediff from genometools gives me some scores that I can use to generate nice clusters, but I am not sure whether that tools really works. The fasta-files that I use contain multiple sequences.

I ran some tests and it looks ok.

1_reference.fna:
>1
TAAGTTACT
>2
TAAGTTACA

2_eq_to_ref.fna:
>1
TAAGTTACT
>2
TAAGTTACA

3_tags_diff.fna:
>1asdfadf
TAAGTTACT
>2asdfasdfffa
TAAGTTACA

4_orde_diff.fna:
>2
TAAGTTACA
>1
TAAGTTACT

5_add_subse.fna:
>1
TAAGTTACT
>2
TAAGTTACA
>2
TTACA

6_point_mut.fna:
>1
AAAGTTACT
>2
TAAGTTACA

7_different.fna:
>1
TAAGTTACT
ATTACCTAA
>2
AAAAAAAAA

Then:

gt genomediff --indexname test *fna
7
1_reference.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
2_eq_to_ref.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
3_tags_diff.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
4_orde_diff.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
5_add_subse.fna 0.206969        0.206969        0.206969        0.206969        0.000000        0.212596        1.349206
6_point_mut.fna 0.199527        0.199527        0.199527        0.199527        0.212596        0.000000        0.569180
7_different.fna 0.794496        0.794496        0.794496        0.794496        1.349206        0.569180        0.000000

I am running the clustering on the difference matrix calculated with genomediff: "These distances are Jukes-Cantor corrected divergence between the pairs of genomes, that is, the number of mutations per base between them."

Currently, we are studying S aureus. We genomes are assembled genomes (three different methods). My guess is that the sequences from the plasmids are present. Furthermore, we do have drug-resistance measured in culture. So, we will be able to compare the genomes and the resistances.

... stuff of original post deleted.
On second thoughts what you might be doing is a template assembly of your genomes. It is a possible interpretation of the 1.2. fasta sequences above (i.e. 1 is template).

A microbial genomics professional would advise also performing a de novo assembly, particularly if you are interested in presence and absence. The reason is that if a gene is present in your query and absent in your template, you will miss it. Again it all depends on what bacteria are being assessed, some are more prone to "genetic islands" than others.

You need a collaborator beyond that.

For the onlookers here bacterial genetic behaviour is very different to eukaryotes and what they get up, switching DNA etc.., to would appear bizzare from a higher eukaryotic world.

You mentioned a comparative analysis of 500 strains: I have worked on this in MRSA. Note, the bacteria is important considering the approach you will adopt.

Anyway, you want a single aligned file to produce a phylogeny, Bayes or likelihood. This is a model of point mutations. Bootstrapping doesn't really help because of SNP differences between isolates.

There is a complex problem. Generic phylogeny can hit the buffers in my opinion because often:

low numbers of SNP differences between isolates.
The error across the genome is poorly defined
Multi-clonal infection ignored and your isolates not cloned.
The other issue is the QC of the genome, I've been fairly astonished by the variation (it affects the tree)

Put it all together and you can get a nice tree but the topology of a given MLST is in my opinion unlikely to be correct. A referee may not bother with this, but they might.

You then map the phenotype against the phylogeny (the tips of the tree) and look for clusters. HOWEVER, there are problems:

Drug-resistant e.g. mec cassette (methicillin-resistant ) are often on plasmids and these are lost during isolation. Albeit they can integrate onto the genome. So an isolate can be drug resistance but you fail to find the gene.
Drug-resistance transmission is unlikely to cluster, its too quick, so you don't see nice tight clusters against a theoretically perfect tree
The best approach is drug-resistance in culture.

Hierarchical clustering, has many meanings. It is used in presence/absence data, I demonstrated the method failed (on bacteria) based on cluster analysis. It is currently revised and used in unsupervised deep (or machine) learning prior a training method.  The question then is what is what are you modelling? Presence/absence (of genes), point mutations, epidemiological data, drug-resistance/non-drug resistance?

My assessment is that given the potential complexity of these biological scenarios and the quality of the question you should seek formal collaboration and this is before considering what bacteria are being assessed. Some of epidemic drug-resistent bacteria rip up this rule book.


How to query a network for a bacterium, specifically Streptomyces caatingaensis or Streptomyces thioluteus? I have a list of gene/protein IDs and I could not query a network using String, BioGrid, Intact and IIS, for Streptomyces thioluteus.

The bacteria I am working with is the Streptomyces caatingaensis. Since it is newer than others, I thought I would find an interaction network looking for Streptomyces thioluteus (which shares one pathway I am analyzing in common, others do not share). As a result, I have this list of proteins, but I could not query any network from the IDs. Sometimes I could not find the bacterium in the tools and when it was found, as in IIS the result was zero interactions.

Is there anything I could use to have some network that could be used as insight for the analysis?

I realized that some tools work better with Gene Names instead of some type of ID.

Using IIS I was capable of creating a network using this list of proteins (Streptomyces caatingaensis) by querying proteins by Genes Names having as base the E. Coli.

The tool has 3 main tabs: proteins, annotation, and interactome.
Search the Gene Names in proteins section, go to Interactome, set the parameters and create the network.

P.S. If you are not used to the tool you should read the research article available on the home page before using the network for research purposes.


How can I extract only insertions from a VCF file? I'm looking to subset a standard VCF file to generate one which only includes insertions (i.e. not indels).

I can get part of the way there with:

bcftools view -v indels <vcf> | awk '{if(length($4) == 1) print}'

However this wouldn't catch an insertion that was part of a multi-allelic record with an insertion and deletion, where the reference length would be greater than 1 bp. Potentially then one way to go is a chain of decomposition, normalisation and then the same reference length filtering — surely there's a better way in one of the many VCF manipulation utils?

A one-liner:

zcat my.vcf.gz |
  perl -ane '$x=0;for $y (split(",",$F[4])){$x=1 if length($y)>length($F[3])}print if /^#/||$x'

or equivalently

zcat my.vcf.gz |
  perl -ane '$x=0;map{$x=1 if length>length($F[3])}split(",",$F[4]);print if /^#/||$x'

For simple VCF operations, I generally recommend to write a script. This may be faster than those using heavy libraries. With a script, you only parse fields you care about; most libraries unnecessarily parse every field.

On a related note, I recommend not to decompose multi-allelic sites unless necessary. Decomposing is tricky, makes VCF harder to parse and to understand and may be a potential source of errors. Here is an example:

#CHROM POS ID REF  ALT               QUAL FILTER  INFO  FORMAT  S1  S2  S3  S4
11     101 .  GCGT G,GCGA,GTGA,CCGT  199  PASS    .     GT      0/1 1/2 2/3 2/4

vt decompose+normalize produces the following VCF:

#CHROM POS ID REF  ALT QUAL FILTER INFO FORMAT  S1   S2   S3   S4
11     101 .  GCGT G   199  PASS   .    GT      0/1  1/.  ./.  ./.
11     101 .  G    C   199  PASS   .    GT      0/.  ./.  ./.  ./1
11     102 .  CGT  TGA 199  PASS   .    GT      0/.  ./.  ./1  ./.
11     104 .  T    A   199  PASS   .    GT      0/.  ./1  1/.  1/.

In theory, you can reconstruct the original VCF from this output. However, it is very challenging for a program to do that. When you compute allele frequency line-by-line, this VCF will give you wrong results. bcftools norm -m- replaces "." with "0". You can get a correct ALT allele frequency from the bcftools output, but a wrong REF allele frequency. Furthermore, vt is also imperfect in that "CGT=>TGA" is not decomposed.

My preferred output is:

#CHROM POS    ID     REF    ALT    QUAL   FILTER INFO   FORMAT S1     S2     S3     S4
11     101    .      GCGT   G,<M>  0      .      .      GT     0/1    1/2    2/2    2/2
11     101    .      G      C,<M>  0      .      .      GT     0/2    2/0    0/0    0/1
11     102    .      C      T,<M>  0      .      .      GT     0/2    2/0    0/1    0/0
11     104    .      T      A,<M>  0      .      .      GT     0/2    2/1    1/1    1/0

Here we use a symbolic allele <M> to represent "another ALT allele". You can calculate the allele frequency by looking at one line, and won't confuse other ALT alleles with REF. bgt can produce such a VCF indirectly. However, it discards all INFO, so is not a practical solution, either.

In summary, it is very difficult to decompose multi-allelic sites. When you get decomposing wrong, your downstream analyses may be inaccurate. Decomposition should be used with caution.


NCBI blast for exact match of a short sequence I'm trying to Blast for exact matches to the sequence: 'ATTGNNNNGCAAACCA' in the human transcriptome using NCBI Blast on its 'refseq_rna' database. However, when I do a basic query I get "No significant similarity found."
Note the N's in the middle of the sequence, if I get rid of those N's, yielding the sequence 'AGCGGATTGCAAAGCAAACCA', I get a match in the 3' UTR of the human MeCP2 gene (which is correct). I don't understand why adding the N's here make it not work.
I looked at the help section and I believe this may be due to shorter sequences being less "statistically significant". Statistical significance of output is completely irrelevant to me, so I set 'Expect Value' to 100000000000000. However, I still get the same result.
I also tried to modify the HTML of the page to allow me to submit a word size of 4 to the form. This went through, but I still got "No significant similarity found".
Could someone help me out with this? I feel like searching for a short sequence shouldn't be this hard.

UPDATE: while procrastinating, I worked up the below throwaway code into a more generally applicable (but still bare-bones) script.

ORIGINAL:
I feel like you might want to use something other than BLAST to do this.
BLAST is not really designed to solve this problem, I think it's just easy to use. Also, your N-containing sequence (16nt) is shorter than the version without Ns (21 nt), which can't be correct if you care about "exact" matches.
As a side note, the idea of exact matches is a bit tricky with Ns involved. The whole point of ambiguity codes is that things won't be exact.
And all of this is on top of the short-string problem- most alignment programs won't even return results that don't meet some minimum alignment score, which is exacerbated a bit by the Ns. For how BLAST handles ambiguity codes, you can see section 7.2 here.
Why not use regular expressions instead to search your sequence against the transcriptome?
# python
import re
from Bio import SeqIO

fwd = r'ATTG[ACGT]{4}GCAAACCA'
rev = r'TGGTTTGC[ACGT]{4}CAAT'  # reverse complement

# reference file
seqfile = SeqIO.parse("your_reference_file.fa", "fasta")  # a transcriptome reference fasta

for seq in seqfile:
    # maybe use re.match() to get coords but easier to show it with findall
    a = re.findall(fwd, str(seq.seq))
    b = re.findall(rev, str(seq.seq))
    if len(a) > 0:
        print(seq.id, a)
    if len(b) > 0:
        print(seq.id, b)

When I put this in a file and run it against a test case with ATTGAAAAGCAAACCA, it works:
$ python exact_match.py
the_contig ['ATTGAAAAGCAAACCA']

Alternately, you could just blast the 16 non-N versions of that sequence and combine them, as you've shown.


Advantages of paired-end sequencing compared to single end One of the advantages of paired end sequencing over single end is that it doubles the amount of data. Another supposed advantage is that it leads to more accurate reads because if say Read 1 (see picture below) maps to two different regions of the genome, Read 2 can be used to help determine which one of the two regions makes more sense. However, I don't understand how this second advantage pertains only to paired end sequencing. With single-end sequencing, couldn't you just use the read that comes right after read 1 to essentially get a longer net read?

In paired-end sequencing, you have an expected fragment length. So if Read 1 maps to two different places, whereas Read 2 only maps to one place. Then you can determine which one is more likely for Read 1 to map to. It cannot be too far away from Read 2 given the expected fragment length.
It is not necessarily true that paired-end sequencing gets twice the amount of data. The total amount of data can be determined when running sequencing. If single-end has 20M reads, I can just stop sequencing at 20M as well for paired-end.


Find covariation of different DNA-binding protein binding between two conditions I have ChIP-seq data of RNAPII and two other factors which we think follow RNAPII during transcription in two different conditions, and I'd like to show that the genes that lose RNAPII also lose the other two factors, and same for gains.
Do you know any methods that can calculate some sort of covariation coefficients?
I know there are methods that calculate covariation of genes between two conditions, but none for different measures on the same gene.

Answer from @stupidwolf converted from comment:
I would try something like this, for every gene you have a logFC of RNAPII in control versus treatment. Do this for other 2 proteins, and you correlate the logFC. In very rough terms, correlation is the scaled form of covariance, so this will work.
If you have 1000 genes, do correlation between 1000 logFC from RNAPII against 1000 logFC from other proteins.
Answer from @m converted from comment:
You use a covariation stat and then perform Pearson's correlation to standardise the response between -1 to 1. An answer of >0.7 is considered robust and >0.8 highly robust. ITs an interesting question and you are right covariation is the way to answerr it


Pulling out a certain gene in a volcano plot I have generated a volcano plot with a differential expression file.

Code for inputing file:

macrophage_list <- read.table("differential_expression_macrophage.csv", header = T, sep = ",")

library(tidyr)
final_df <- df %>% pivot_longer(., -c(Feature.ID, Feature.Name), names_to = c("set",".value"), names_pattern = "(.+)_(.+)")

# A tibble: 80 x 6
   Feature.ID Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
   <fct>      <fct>        <chr>           <dbl>            <dbl>            <dbl>
 1 a          A            Cluster.1    0.000961            0.292           1
 2 a          A            Cluster.2    0.000902            0.793           1
 3 a          A            Cluster.3    0.00181             1.46            0.758
 4 a          A            Cluster.4    0.000642            0.269           1
 5 b          B            Cluster.1    0.000320            1.95            0.910
 6 b          B            Cluster.2    0.00180             4.77            0.154
 7 b          B            Cluster.3    0                   2.19            1
 8 b          B            Cluster.4    0                   1.66            1
 9 c          C            Cluster.1    0.00128            -2.01            0.0467
10 c          C            Cluster.2    0.00632             0.352           1
# … with 70 more rows

output of head(final_tumor)

# A tibble: 6 x 6
  Feature.ID        Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
  <fct>             <fct>        <chr>           <dbl>            <dbl>            <dbl>
1 ENSG00000227232.5 WASH7P       Cluster.1     0                   1.50            1
2 ENSG00000227232.5 WASH7P       Cluster.2     0                   1.73            1
3 ENSG00000227232.5 WASH7P       Cluster.3     0                   1.77            1
4 ENSG00000227232.5 WASH7P       Cluster.4     0.00114             4.30            0.293
5 ENSG00000227232.5 WASH7P       Cluster.5     0                   2.15            1
6 ENSG00000227232.5 WASH7P       Cluster.6     0                   1.22            1

output of tail(final_tumor)

# A tibble: 6 x 6
  Feature.ID        Feature.Name set        Mean.Counts Log2.fold.change Adjusted.p.value
  <fct>             <fct>        <chr>            <dbl>            <dbl>            <dbl>
1 ENSG00000210196.2 MT-TP        Cluster.6       0.0699          -0.202           0.790
2 ENSG00000210196.2 MT-TP        Cluster.7       0.0801           0.0386          1
3 ENSG00000210196.2 MT-TP        Cluster.8       0.0711           0.0875          1
4 ENSG00000210196.2 MT-TP        Cluster.9       0.0152          -2.31            0.00127
5 ENSG00000210196.2 MT-TP        Cluster.10      0.0147          -2.30            0.00612
6 ENSG00000210196.2 MT-TP        Cluster.11      0.122            0.762           1

Code for generating volcano plot:

library(ggplot2)
library(ggrepel)
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05),
                  aes(label = Feature.Name))

Now, I want to pull out a certain gene, Casp14, from the list and box it on the plot. How do I do that?

As you did for labeling genes with an adjusted p value below 0.05, you can subset your dataset for keeping only rows corresponding to "Casp14":

library(ggplot2)
library(ggrepel)
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05),
                  aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"),
                  aes(label = Feature.Name), color = "red")

With this code, you should see now the labeling of the gene of interest (Casp14) in red on your volcano plot.


well identifiers in a 1536-well mtp plate There's a commonly accepted convention for naming the wells in a 384-well-plate:

A01..A24
B01..B24
..
P01..P24

It's concise and unambiguous and thus very useful.
For 1536-well plates (32x48) it's more difficult because there's more rows than letters in the alphabet. It seems there's no standard.
Which nomenclatures are out in the field and which is the best?

and here's a scheme I made up myself:

AA01..AA48
..
AZ01..AZ48
BA01..BA48
..
BF01..BF48

Pros:

alphabetic sorting of well ids does what you expect
it's clear what's the row and what's the column
it's easy to distinguish well ids of 384-well plate from those of a 1536-well plate

Cons:

it's not used anywhere AFAIK


bedtools feature out of bounds I'm trying to count the GC in regions leading up to polyA clusters. I found the bedtools command nuc, one of the outputs of which is the GC content.

I took the mouse polyA cluster database from here, and the mm10 primary assembly soft masked from here. I imported the polyA clusters bed file to R. I then got a list of all the polyA clusters in the genes I'm interested in. For positive strands on polyA I took the start position, took 100 away from it and made that the start of the region of interest, then made the start of the polyA signal the end. For things on the negative strand I added 100 to the end of the polyA signal. I then exported the results to a bed file format. The code I used to do this is below.

When I try to run:

bedtools nuc -fi Mus_musculus.GRCm38.dna_sm.primary_assembly.fa -bed PolyA.bed

Every signal returns some variant of: Feature (chr10:21,080,315-21,080,415) beyond the length of chr10 size (0 bp).  Skipping.

I can however look at several of these regions in the UCSC genome browser without a problem, at least I don't think I had a problem.

Can anyone point me to what's going wrong or to another way I can solve this.

GetRanges <- function(FilteredClusters){
  for (x in 1:nrow(FilteredClusters)){
    if(FilteredClusters$Strand[x] == '+'){
  FilteredClusters$RegionEnd[x] <- FilteredClusters$ClusterStart[x]
  FilteredClusters$RegionStart[x] <- FilteredClusters$RegionEnd[x] - 100
}
else if(FilteredClusters$Strand[x] == '-')
    {
      FilteredClusters$RegionStart[x] <- FilteredClusters$ClusterEnd[x]
      FilteredClusters$RegionEnd[x] <- FilteredClusters$RegionStart[x] + 100
    }
  }
    ClusterCut <- select(FilteredClusters, 'Chr', 'RegionStart', 'RegionEnd')

    return(ClusterCut)
}

There is no chromosome named chr10, but there is one named 10. The problem is that you're mixing chromosome name systems with a tool that, like most, isn't able to handle that. Remove all instances of chr in your BED file and fix things like chrM -> MT.


How can I apply proportional (p) distances (Nucleotide) using bioPython I am working on phylogenetic tree generation. I used BioPython and ClustalW 2.x.x for this purpose. I have generate the tree using BioPython but when I try to generate tree using "MegaSoftware GUI" tool my tree does not matched with the output of MegaSoftware.

When I use "No. of difference" model then results are matched with biopython tree Like below image shown.

But whenever I choose "p-distance" model my results are not matched.

Output MegaSoftware using p-distance

Output BioPython

Here is my code for generating phylogenetic tree:

def conversion_fasta_dnd(clustalw, filepath, type, clustering, output, pwmatrix):
    """
        ClustalW tool to conversion the file format.
        :param clustalw: path of the executable programm.
        :param filepath: path of `fasta` extension file.
        :param type: dna, protein
        :param clustering: upgma, nj
        :param output: clustal, phylip, fasta
        :param pwmatrix: blosum, pam, id
        :return: dnd, aln (files)
        """
    # esblish the full path to the clustalW program
    clustalw_exe = clustalw

    # check to make sure the program is there
    assert os.path.isfile(clustalw_exe), "Clustal W executable missing"

    # create the appropriate command line
    clustalw_cline = ClustalwCommandline(clustalw_exe,
                                         infile=filepath,
                                         type=type,
                                         clustering=clustering,
                                         output=output,
                                         pwmatrix=pwmatrix)

    # print the command line
    print(clustalw_cline)

    # execute the command
    stdout, stderr = clustalw_cline()

# This function draw the tree on image file
def simple_tree(filepath):
    """
    Generate the simple tree image.
    :param filepath: path of dnd file format.
    """
    def get_label(leaf):
        return leaf.name

    file = filepath
    tree = Phylo.read(file, 'newick')
    Phylo.draw(tree, label_func=get_label, do_show=False, branch_labels=lambda c: c.branch_length)

    # axis management
    pylab.axis('on')
    pylab.savefig(OUTPUT_IMAGE, format='png', bbox_inches='tight', dpi=100, bbox_extra_artists=[])

In above code I have used two things first BioPhython lib and the other is Clustal W for align the sequences.

Could you please clarify what p-distances are?
I only have these reference for p-distance:

p-distance (Nucleotide).
Proportional (p) distances

I need to add p-distance model in my code, someone please help me out How can I do that.

Based on the sites you mention, it seems that p-distances are pairwise distances between aligned sequences obtained by counting the number of differences and dividing by the length of the alignment. This is the proportion of differences between two sequences.

I'm not sure exactly what you mean by "add p-distance model in my code", but here is how you could compute those distances from an existing alignment (here I use as an example an alignment I have in fasta format):

from Bio import AlignIO

def p_dist(seq1, seq2):
    ali_len = len(seq1)
    assert ali_len == len(seq2)
    return sum(0 if nuc1 == nuc2 else 1 for (nuc1, nuc2) in zip(seq1, seq2)) / ali_len

ali = AlignIO.read("most_targeted_histones/transcript_alignments/H4_top4_loss_spliced_NT.fa", "fasta")
nb_seq = len(ali)
distances = {}
for i in range(nb_seq - 1):
    rec1 = ali[i]
    j = i + 1
    while j < nb_seq:
        rec2 = ali[j]
        distances[(rec1.id, rec2.id)] = p_dist(rec1.seq, rec2.seq)
        j += 1

print(*distances.items(), sep="\n")

Output:

(('ZK131.4_spliced', 'ZK131.8_spliced'), 0.0)
(('ZK131.4_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.4_spliced', 'ZK131.1_spliced'), 0.0)
(('ZK131.8_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.8_spliced', 'ZK131.1_spliced'), 0.0)
(('T10C6.14_spliced', 'ZK131.1_spliced'), 0.10897435897435898)

If you want to build a tree from those distances, I'm not sure you can do this with BioPython. It seems however that there are some possibilities with DendroPy.

You would need to write the distances in a suitable format, and to choose the method between UPGMA and neighbour-joining.


Error When Using biocLite as an installer in rpy2 python library I am running into an error when I tried to run the following lines of code, the error at line 3:
1 from rpy2.robjects.packages import importr
2 base = importr('base')
3 base.source("http://www.bioconductor.org/biocLite.R")
4 biocinstaller = importr("BiocInstaller")
5 biocinstaller.biocLite("seqLogo")

rpy2.rinterface_lib.embedded.RRuntimeError: Error: With R version 3.5 or greater, install Bioconductor packages using BiocManager; see https://bioconductor.org/install
How can I use BiocManager installer instead of biocLite  and reflect that on my code? I couldn't find any documentation or source that talk about this error from within Python.
Thank you in advance

You're using an outdated version of R. You can use the R 4.0 from conda-forge channel, just create an environment with r-base and then install rpy2, the latest version compatible with R 4.0 is available through pip
conda create -n rpy2 -c conda-forge r-base
conda activate rpy2
pip install rpy2

Now you can use rpy2 with R 4.0 in your python code.
import rpy2.robjects as robjects
robjects.r('install.packages("BiocManager", '
           'repos="http://cran.r-project.org")')
robjects.r('BiocManager::install("seqLogo")')
robjects.r('library(seqLogo)')

Remember deactivate your env with conda deactivate.


How to design DESeq2 LRT model with individuals nested in 2 levels? We have a complicated experimental design that we would like to perform LRT analysis for. Our main goal is to discover significant genes for the "Injection:Social" interaction term across the entire dataset by removing it from the LRT reduced model, and as a bonus we are also interested in discovering significant genes for that interaction term for each respective brain region.
Sample  Injection   Social  Region  Individual  ind.n
HY06    L   ISO HY  S06 S1
NST6    L   ISO NS  S06 S1
TN06    L   ISO TN  S06 S1
HY08    L   ISO HY  S08 S2
NST8    L   ISO NS  S08 S2
TN08    L   ISO TN  S08 S2
HY30    L   KF  HY  S30 S1
NST30   L   KF  NS  S30 S1
TN30    L   KF  TN  S30 S1
HY32    L   KF  HY  S32 S2
NST32   L   KF  NS  S32 S2
TN32    L   KF  TN  S32 S2
HY64    L   KFC HY  S64 S1
NST64   L   KFC NS  S64 S1
TN64    L   KFC TN  S64 S1
HY65    L   KFC HY  S65 S2
NST65   L   KFC NS  S65 S2
TN65    L   KFC TN  S65 S2
HY19    L   NF  HY  S19 S1
NST19   L   NF  NS  S19 S1
TN19    L   NF  TN  S19 S1
HY24    L   NF  HY  S24 S2
NST24   L   NF  NS  S24 S2
TN24    L   NF  TN  S24 S2
HY05    S   ISO HY  S05 S1
NST5    S   ISO NS  S05 S1
TN05    S   ISO TN  S05 S1
HY12    S   ISO HY  S12 S2
NST12   S   ISO NS  S12 S2
TN12    S   ISO TN  S12 S2
HY31    S   KF  HY  S31 S1
NST31   S   KF  NS  S31 S1
TN31    S   KF  TN  S31 S1
HY34    S   KF  HY  S34 S2
NST34   S   KF  NS  S34 S2
TN34    S   KF  TN  S34 S2
HY62    S   KFC HY  S62 S1
NST62   S   KFC NS  S62 S1
TN62    S   KFC TN  S62 S1
HY63    S   KFC HY  S63 S2
NST63   S   KFC NS  S63 S2
TN63    S   KFC TN  S63 S2
HY04    S   NF  HY  S04 S1
NST4    S   NF  NS  S04 S1
TN04    S   NF  TN  S04 S1
HY20    S   NF  HY  S20 S2
NST20   S   NF  NS  S20 S2
TN20    S   NF  TN  S20 S2

My first attempt was building simple full (m1) and reduced (m2) models that gets directly at our question of interest but doesn't control for nested individuals.
m1 <- model.matrix(~ Region + Social * Injection, colData_filt)
m2 <- model.matrix(~ Region + Social + Injection, colData_filt)

We want to control for individual/batch effects, which is nested within both "Injection" and "Social" but not region, as we have three brain regions per individual. I followed the example in the DESeq2 manual for creating a term (ind.n) distiguishing individuals nested within groups, but now I'm not sure how to create the full and reduced model given that I have one more level than the example.
I've tried a really elaborate full model (m1) with the interaction term of interest (Injection:Social) removed for the reduced model (m2), but I'm not sure this is correct based on our design.
m1 <- model.matrix(~ Injection + Injection:ind.n + Injection:Social + Injection:Region + Social + Social:ind.n + Social:Region + Region, colData_filt)
m2 <- model.matrix(~ Injection + Injection:ind.n + Injection:Region + Social + Social:ind.n + Social:Region + Region, colData_filt)

I'm assuming this is wrong, but even if this was by some miracle the correct formulation, would there be a way to extract genes that explain the "Injection:Social" interaction term for separate brain regions?
As a work-around, I subsetted the data by region and ran three separate LRT analyses for each subset and compared the results. While this simplified the model to look like the first example above, I worry that we lose some power by ignoring the fact we have multiple brain region samples from single individuals across the dataset.
Any guidance is much appreciated. Thanks in advance

From what I can gather, you want to account for the effect of individual, nested within region. That is, you want to see after accounting for these, is there a consistent effect for Injection:Social across all conditions.
So you set up the model like this:
m1 <- model.matrix(~ ind.n*Region + Injection + Social + Injection:Social,data=..)

The last term should be Injection:Region and you can just use the waldTest (default) in DESeq2 for this term.
What does the terms do? ind.n*Region is the equivalent of ind.n + Region + ind.n:Region , and with this you effectively get an effect for every region in every individual.
Why don't we need the Injection:ind.n or Social:ind.n or Social:Region. These terms indicate the effect of Injection or Social can vary by individuals or regions. Most likely introducing too many parameters when you are interested in a common effect. Also you do not have the replicates or samples to distinguish this effect from region or other effects.
Since you provided an example we can run DESeq2 and you can see how the results look like:
mat = counts(makeExampleDESeqDataSet(n=1000,m=48))
dds = DESeqDataSetFromMatrix(mat,df,~ ind.n*Region + Injection + Social + Injection:Social)

dds = DESeq(dds)

resultsNames(dds)
 [1] "Intercept"            "ind.n_S2_vs_S1"       "Region_NS_vs_HY"
 [4] "Region_TN_vs_HY"      "Injection_S_vs_L"     "Social_KF_vs_ISO"
 [7] "Social_KFC_vs_ISO"    "Social_NF_vs_ISO"     "ind.nS2.RegionNS"
[10] "ind.nS2.RegionTN"     "InjectionS.SocialKF"  "InjectionS.SocialKFC"
[13] "InjectionS.SocialNF"

The terms you need are "InjectionS.SocialKF","InjectionS.SocialKFC", "InjectionS.SocialNF", and you can look at each of them:
head(results(dds,name="InjectionS.SocialNF"))
log2 fold change (MLE): InjectionS.SocialNF
Wald test p-value: InjectionS.SocialNF
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat
             <numeric>          <numeric>         <numeric>          <numeric>
gene1  9.9811166787259   1.25304112986447 0.819806376919295   1.52845984752303
gene2 30.3449455820337 0.0329442893152027 0.705199688255367 0.0467162562092241
gene3 3.83223545055379    1.0281136369045  1.64095596190233  0.626533350543196
gene4  11.232305747171  0.595738624408923   0.8243883031544  0.722643227868976
gene5 6.70950627004097  0.756449993378065   1.0631622863378  0.711509430967263
gene6 26.1431134888287 -0.854784518963918 0.625714541243558  -1.36609342219393
                 pvalue              padj
              <numeric>         <numeric>
gene1 0.126398405431826 0.978671002658464
gene2 0.962739373909937 0.999897888026606
gene3 0.530965168963838 0.978671002658464
gene4 0.469899103018657 0.978671002658464
gene5 0.476768608734069 0.978671002658464
gene6 0.171909642630577 0.978671002658464

As mentioned, you can do a LRT if you want to test all the Injection:Social term interactions terms at one go, that is, the null hypothesis is that all of them are zero:
dds = nbinomLRT(dds,reduced=~ ind.n*Region + Injection + Social)
results(dds)

Usually the individual terms are more intuitively than testing all of them are zero, but you might have a special need for this.


How to calculate overall reference coverage with MUMmer? Is the MUMmer suite capable of calculating reference sequence coverage statistics for all query sequences collectively? It would be possible to achieve by parsing the output of nucmer / show-coords / show-tiling but it seems like there should be a better way. I currently do this using a sensitive read mapper, samtools depth and some scripting.

To clarify, I'd like to know the reference coverage achieved using all of the query sequences (i.e. the whole de novo assembly, in my case).

I believe all you need to do is to run dnadiff from MUMmer. That will run a comparison and output a number of useful metrics.


How to convert species names into common names? I’m trying to find common names from a list of scientific names (not all will have them though).

I was attempting to use taxize in R but it aborts if it doesn’t find an entry in EOL and I don’t know a way around this other than manually editing the list - in which case I might as well check all this way.

Does anyone have a better way of doing this?

use the NCBI taxon dump under ftp://ftp.ncbi.nih.gov/pub/taxonomy

in taxdmp.zip you'll find all the names for a given NCBI taxon

$ grep -w ^9606 names.dmp
9606    |   Homo sapiens    |       |   scientific name |
9606    |   Homo sapiens Linnaeus, 1758 |       |   authority   |
9606    |   human   |       |   genbank common name |
9606    |   man |       |   common name |


Core genome too short I'm trying to create a snp distance matrix from M tuberculosis isolates to try and infer transmission networks. My plan was to use snippy to make a core genome alignment and then snp-dists to create the matrix but have hit a snag with creating a core genome alignment. When I try it with all 972 of my samples, snippy halts when trying to run snp-sites, saying Warning: No SNPs were detected so there is nothing to output. The issue is that with all the samples included, the core genome size is 0.
When looking at the alignment statistics provided by snippy, it looks like some isolates have very low numbers of aligned base pairs (??due to low coverage or contamination). I have progressively filtered out samples with low numbers of aligned base pairs and can then get snippy/snp-sites to work. As I increase my filter stringency, the resulting core genome is still short. Eg Filtering out those with aligned base pairs of <90% of the reference had a core genome of only 5747bp (from a reference of 4.4 million bp), while excluding 88 (9%) of my samples.
snippy's author recommends using it's output file core.txt to figure out which samples are the "bad" outliers. That file provides 1) the length of the reference, 2) the number of aligned base pairs, 3) the number of unaligned base pairs, 3) the number of variant sites, 4) the number of heterogenous sites, 5) the number of masked sites, and 6) the number of low coverage sites.
#>   ID                        LENGTH ALIGNED UNALIGNED VARIANT   HET MASKED LOWCOV
#>   <chr>                      <dbl>   <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>
#> 1 R15795_CATCAAGT_S34_L006 4411532 4192432      6818     755   428 209178   2676
#> 2 R15842_GTCTGTCA_S49_L006 4411532 4187344     10304     784   441 209178   4265
#> 3 R15876_CGCTGATC_S36_L006 4411532 4176662     18992    1292   660 209178   6040
#> 4 R15951_ATTGGCTC_S7_L002  4411532 4170649     14732    1281   980 209178  15993
#> 5 R16019_TGGAACAA_S78_L001 4411532 4190733      8132     712   715 209178   2774
#> 6 R16046_GACTAGTA_S6_L002  4411532 4186069      4309    1257  1084 209178  10892

My question is: what heuristic would you use to filter out samples prior to creating a core genome?
And secondarily, what would a reasonable core genome size be for Mtb?

Targeting Illumina reads with >70X depth-of-coverage is a good place to start.  Once a robust database is established lower coverage isolates can be used effectively too.
More specific to individual SNPs... filtering on VCF QUAL value, Map Quality and/or AC value have shown to be useful.
Core genome will likely be around 1,000 SNPs, but it will depend on the lineage the reference belongs to and the aligned sample lineage.  If both are from the same lineage you're likely to be less than 700 SNPs.  When from different lineages SNP counts will obviously be higher in addition to the alignment error seen.
Pertaining to filtering individual SNPs...
vSNP is a tool developed for SNP analysis.  It can be applied to different specie types but was originally developed for use with TB complex isolates.  Code detail is here.
vSNP is a 2-step process.  First, generating VCF files.  Second, using those VCF files to output SNP tables and trees.  The goal of the script is to achieve a high resolution genotyping analysis that can be validated with minimal effort.  Some steps are included to help achieve this such as reference selection, automatic grouping and position filtering.
Pseudocode:
For set of VCF files created from the same reference...
def find_positions()
    AC=2
    QUAL > 300

For each found position...
def get_snps()
    AC=2 (solid call)
    QUAL > 50
        THEN --> SNP
    AC=1 (mix call)
    QUAL > 150
        THEN --> AMBIGIOUS
    AC=2
    QUAL >= 50
       THEN --> REF
    QUAL < 50
       THEN --> N
    ALT == "None"
       THEN --> -

First find all possible positions that could be informative in all samples/VCF files: find_positions(), then define each position to be included in the alignment: get_snps().
Then it can be convenient to filter consistently poor alignment regions such as PPE/PGRS.


Writing a perl script to holding information for two genes Basically I have a perl script in which I have an array (where each element of the array references a hash) and need to be able to print the array with a dumper function.
Thus I need to be able to split the $line on white spaces and save into variables and then construct anonymous hash and push into @genes:

#!/usr/bin/perl
use strict;
use warnings;
use Data::Dumper; # print complex data structure
my @genes; # declare the array
while(<>) { # this means that as long as lines come from the pipe we keep going
  my $line = $_; # a line that come from the pipe (we go line by line)
  next unless $line =~ /^\d+/; # skip lines except those reporting genes
  <split the $line on white spaces and save into variables>
  <anonymous hash and push into @genes>
}
print Dumper(\@genes);
exit;

I don't  even know how to approach this, this is my first time perl scripting and I am really confused. Any insight would be extremely helpful.

So far I have:

my $gene = {'id' => $id, 'start' => $start, 'end' => $end, 'frame' = $frame, 'score' => $score};
my @line_split = split(' ', $line);
my $id = $line_split [0], $start = $line_split [1], $end = $line_split [2], $frame = $line_split [3], $score = $line_split [4]

Perhaps this small test script will help demonstrate some of the principles:

#!/usr/bin/env perl

use strict;
use warnings;
use Data::Dumper;

my @arr;

# define the current line
my $line = "foo\tbar\tbaz\n";

# remove the trailing newline character
chomp $line;

# put the line's elements into an anonymous array
my ($a, $b, $c) = split('\t', $line);

# insert key-value pairs into anonymous hash
my $hash = {a => $a, b => $b, c => $c};

# add the hash table to the array
push @arr, $hash;

# print the array
print Dumper \@arr;

Run this script and see what it prints out.


Missing genes and normalisation of RSEM output using EBSeq Without going into too much background, I just joined up with a lab as a bioinformatics intern while I'm completing my masters degree in the field. The lab has data from an RNA-seq they outsourced, but the only problem is that the only data they have is preprocessed from the company that did the sequencing: filtering the reads, aligning them, and putting the aligned reads through RSEM. I currently have output from RSEM for each of the four samples consisting of: gene id, transcript id(s), length, expected count, and FPKM. I am attempting to get the FASTQ files from the sequencing, but for now, this is what I have, and I'm trying to get something out of it if possible.

I found this article that talks about how expected read counts can be better than raw read counts when analyzing differential expression using EBSeq; it's just one guy's opinion, and it's from 2014, so it may be wrong or outdated, but I thought I'd give it a try since I have the expected counts.

However, I have just a couple of questions about running EBSeq that I can't find the answers to:

1: In the output RSEM files I have, not all genes are represented in each, about 80% of them are, but for the ones that aren't, should I remove them before analysis with EBSeq? It runs when I do, but I'm not sure if it is correct.

2: How do I know which normalization factor to use when running EBSeq? This is more of a conceptual question rather than a technical question.

Thanks!

Yes, that blog post does represent just one guy's opinion (hi!) and it does date all the way back to 2014, which is, like, decades in genomics years. :-) By the way, there is quite a bit of literature discussing the improvements that expected read counts derived from an Expectation Maximization algorithm provide over raw read counts. I'd suggest reading the RSEM papers for a start[1][2].

But your main question is about the mechanics of running RSEM and EBSeq. First, RSEM was written explicitly to be compatible with EBSeq, so I'd be very surprised if it does not work correctly out-of-the-box. Second, EBSeq's MedianNorm function worked very well in my experience for normalizing the library counts. Along those lines, the blog you mentioned above has another post that you may find useful.

But all joking aside, these tools are indeed dated. Alignment-free RNA-Seq tools provide orders-of-magnitude improvements in runtime over the older alignment-based alternatives, with comparable accuracy. Sailfish was the first in a growing list of tools that now includes Salmon and Kallisto. When starting a new analysis from scratch (i.e. if you ever get the original FASTQ files), there's really no good reason not to estimate expression using these much faster tools, followed by a differential expression analysis with DESeq2, edgeR, or sleuth.

1Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN (2010) RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics, 26(4):493–500, doi:10.1093/bioinformatics/btp692.

2Li B, Dewey C (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12:323, doi:10.1186/1471-2105-12-323.


Why does samtools mpileup sometimes include ref bases (other than ',' or '.')? This is my first post here. I can think of no way of giving an easily reproducible example, as per the stack ethos. SO apologies in advance and any feedback on question format appreciated.

I have called samtools mpileup (Version: 1.9 (using htslib 1.9)) on some oxford nanopore data (mapped to GRCh38).

I am seeing some instances like below (grouped with pythons Counter), here reference is T at this position:

 Counter({',': 34, '.': 33, 'c': 9, 'g': 6, 'a': 4, 't': 2, '-:ACTC': 1, '-:ACTCGG': 1, '-:5': 1, '$': 1, '*': 1})

Why is 't': 2' in there?

flags used: -A -B -Q 1 -R -x

Why is 't': 2' in there?

How else do you think the pileup should show you that two reads (in the reverse direction) have a T in that position?