Are there any tools for discovering tandem repeats around a variant in a VCF file? I know for example PolyX from Pierre's Jvarkit to find the number of repeated REF bases (homopolymers) around a position in a VCF file as described here exists.

I would like to know however if there tools that describes tandem repeats around a reported variant (TR) up to 6 repeated bases. I can imagine 3 different scenarios:

Scenario 1:

chrK    posY    C     G

TR of (ATC)
____________v
ATCATCATCATCGATGATAC

Scenario 2:

chrK    posY    A     T

             TR of (AACT)
            v____________
CAGTCGATTGGATAACTAACTAACT

Scenario 3:

chrK    posY    C     A

TR of (AG)    TR of (AG)
____________v__________
AGAGAGAGAGAGCAGAGAGAGAG

If you know any tools in mind that can be used for VCFs, I'd really appreciate citing them.

Thanks in advance.

After researching a bit, I found a couple of tools that take as input BAM files. I did not test them yet:

GangSTR
lobSTR
hipSTR
Tandem Repeat Finder (thanks to Michael G.'s answer)

If I find more tools, I will update this post. After testing them, I will come back to share my feedback.


How to compute the chance of failing to detect a gene given the detection limit of a protocol In Shapiro et al., when discussing about loss of molecules as source of error in single-cell sequencing, it is written that:

  Another source of error is losses, which can be severe. The detection limit of published protocols is $5$–$10$ molecules of mRNA. If, as seems likely, the limit of detection is primarily determined by losses during sample preparation, this would indicate that $80$–$90\%$ of mRNA was lost. Or, to put it the other way around, a $90\%$ loss leads to an approximately $50\%$ chance of failing to detect a gene that is expressed at a level of seven mRNA molecules (from the binomial distribution).

How is this probability computed using the binomial distribution? I thought that $90\%$ loss corresponds to $5$ detected molecules, and I assume that $k=7$ for the binomial calculation, but I am unable to go further.

A 90% loss can be rephrased as a 10% chance of detecting anything. So what we want to find is the probability of detecting 0 molecules, when we start with 7 and have 10% probability of success. Once can do that in R as follows:

> pbinom(0, 7, 0.1)
0.4782969

So ~50%, as they stated. I suspect that part of the confusion arises from the fact that the detection limit is due to lowly-expressed genes/transcript being heavily affected by this loss. So the probability of detecting a single molecule out of 4 original molecules (assuming 90% mRNA loss) is 34%, for 3 molecules it's 27%, for 2 it's 19% and for 1 it's 10%. I think the threshold of 5 is mentioned more because it's a nice round number than there's anything particularly different in detecting a gene with 4 vs. 5 molecules (34 vs. 41% probability).


ipfam or pfam domain-domain interactions I'd like to download domain-domain interactions that were in ipfam (which doesn't seem to be up-and-running anymore).

for an example of what I would like to be able to download via ftp or another source, see this list of interactions with Death domain:
https://pfam.xfam.org/family/PF00531#tabview=tab8

Is this available in one of these files in the current release??
ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release

I think you're looking for ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/database_files/pfamA_interactions.txt.gz. Based on the file's name, that's a list of domain interactions. The file has two columns, so presumably each line represents an interaction.

Here's what I get for the Death domain (PF00531):

$ zgrep PF00531 pfamA_interactions.txt.gz
PF00036 PF00531
PF00531 PF14786
PF00531 PF00531
PF00791 PF00531
PF01335 PF00531
PF13499 PF00531
PF13833 PF00531
PF14786 PF00531

Those are 8 interactions, which matches what the PFam link you gave in your question shows. If you want to convert the accessions to something more informative, you could try Pfam's REST API and do something like:

$ zgrep PF00531 pfamA_interactions.txt.gz |
    while read a b; do
    name1=$(wget -O - "https://pfam.xfam.org/family/${a}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    name2=$(wget -O - "https://pfam.xfam.org/family/${b}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    printf "%s\t%s\t%-20s\t%-25s\n" "$a" "$b" "$name1" "$name2";
done
PF00036 PF00531 EF hand                 Death domain
PF00531 PF14786 Death domain            Tube Death domain
PF00531 PF00531 Death domain            Death domain
PF00791 PF00531 ZU5 domain              Death domain
PF01335 PF00531 Death effector domain   Death domain
PF13499 PF00531 EF-hand domain pair     Death domain
PF13833 PF00531 EF-hand domain pair     Death domain
PF14786 PF00531 Tube Death domain       Death domain

Obviously, if you need anything more fancy than that, it's time to write a proper parser using Perl or Python or any other language that can understand XML.


Theoretical limit of human genome compression How small can a compressed file containing the human genome be?
I'm aware that this question cannot actually be answered, since it is asking for the Kolmogorov complexity of the human genome, which is not computable. So reformulating:
What is the largest known lower bound that we can put on the size of a file containing the human genome?
An example answer could be: on average we know that the human genome contain at least this amount X of randomness, therefore it is not possible to do better than Y bytes.
I'm agnostic on which human genome, but I would say without polymorphisms to make it simpler.
The compressed sequence would need to be stand alone, without using a reference (I asked a question on that topic here: Compressing the human genome to few megabytes).
There should be a program able to compress and decompress the genome in a lossless way. Apart from that, any kind of modelling and coding is allowed.

TLDR; For the most modern CHM13 assembly (v1.1; Complete T2T reconstruction of a human genome; https://github.com/marbl/CHM13), the best compression without any use of a reference would be between 624 MB and 641 MB.
Upper estimates of genome compressibility can be made using standard compressors (see, e.g., this comparison). In my experience, xz can achieve great compression ratios on genomic sequences, especially when used with the maximum compression level (-9) and a single thread (-T1). Importantly, the original FASTA file needs to be cleaned – linebreaks hamper compression as they act as noise; the reformatting can be done by seqtk seq -U (which also converts all nucleotides to the upper case). It the case of the CHM13 assembly, xz compresses the human genome to 641,369,600 bytes.
Lower estimates are trickier. The most straightforward approach is through the number of distinct k-mers. Since k-mer sets contain less information than the original assemblies, we can use them for determining the lower bounds of genome compressibility. The lower bound on the compressibility of k-mer sets is 2 bits per distinct k-mer (under certain assumptions that I believe are satisfied in our case). The CHM13 human genome assembly contains 2,497,988,809 distinct canonical 31-mers, which corresponds to 624,497,203 bytes.
Overall, xz does a great job at compressing single human genomes, and the resulting size is very close to the optimum. The full analysis from this post can be found on http://github.com/karel-brinda/human-genome-compression. The resulting files are here: https://zenodo.org/record/4776580.


*very* unbalanced group sizes for DE I downloaded some publicly available RNA-seq data and want to compare those samples carrying a mutation (~4) against the rest (~800!). I ran both EdgeR and DESeq2, and the first results in an asymmetric volcano plot: skewed in one of the sides, meaning that it results in very large logFC that have large (NS) adjusted p-values. So next, in order to see if this asymmetry is persistant, I ran several random partitions of the data as groups (but keeping them highly unbalanced), obtaining the same behavior. Here an example of what I mean (I run it with 4, 20 and 100 to show how this effect decreases as I increase the small class size)

*I tried EdgeR-robust and I still get the same effect.

I do not see this effect when running DESeq2, only very few DE genes.  What method would you recommend using when comparing such different class sizes?

Random partitioning (as you have done) seems like a reasonable thing to do to work out what's going on. From the results you have suggested, it looks like DESeq2 is performing better than EdgeR, so use DESeq2. Can you show the volcano plots and the partition sizes?

Note: I'm a little bit biased in this, because DESeq2 is what I prefer to use (particularly where Salmon is not appropriate).


State of the art in predicting Translation Initiation Sites I'm working on a university project of predicting Translation Initiation Sites in human DNA. I searched the net for papers and documentation to get guidelines and inspiration, but I feel uncertain that I was able to find the state of the art in solving this problem, so I thought of asking the community here: What's the state of the art method of predicting Translation Initiation Sites?

(I'm also posting the question to help this nascent Q&A site by participating, please close it if it completely misses the mark of being a good question.)

Ribosome Profiling experiments profile the positions of ribosome-protected mRNA. a cycloheximide treatment fixes the translating ribosomes. A naive approach to define Translation Initiation Sites (TIS) would be to look at positions where you see distinct 'peaks' in the distribution of these fragments since the ribosome density will be higher at the start of the CDS (at the TIS).

TISDB is a database comprised of GTI-Seq global mapping of TIS codons at nearly single-nucleotide resolution. GTI-Seq uses two 'fixers' that fix ribosomes where one of them specifically fixes the initiating ribosomes.


How to find genomic overlaps faster? I was trying to find a way to find overlaps between two genomic ranges. I found a post at Biostars but it couldn't consider the chromosome information. For example:

library(IRanges)
df1 = data.frame(chr = c("chr1", "chr12"), start = c(10000, 10000), end = c(20000, 20000))
df2 = data.frame(chr = rep("chr1", 4), posn = c(100, 12000, 15000, 250000), x = rep(1, 4), y = rep(2, 4), z = rep(3, 4))

df1.ir = IRanges(start = df1$start, end = df1$end, names = df1$chr)
df2.ir = IRanges(start = df2$posn, end = df2$posn, names = df2$chr)

df1.it = NCList(df1.ir)

overlap_it = findOverlaps(df2.ir, df1.it)
print(overlap_it)

It returns the overlap from both chromosomes. I am wondering if there is a way to keep the chromosome information intact.

Just wanted to follow up to a question asked at Biostars. My solution and suggestion when dealing with genomic ranges.

library(GenomicRanges)
df1 = data.frame(chr = c("chr1",  "chr12"), start = c(10000,  10000), end = c(20000, 20000))
df2 = data.frame(chr = rep("chr1", 4), posn = c(100, 12000, 15000, 250000), x = rep(1, 4), y = rep(2, 4), z = rep(3, 4))

df1.ir = IRanges(start = df1$start, end = df1$end, names = df1$chr)
df2.ir = IRanges(start = df2$posn, end = df2$posn, names = df2$chr)

df1.it = NCList(df1.ir)

overlap_it = findOverlaps(df2.ir, df1.it)
print(overlap_it)

# define a genomic interval tree (git) for faster search

df1.gr = GRanges (IRanges(start = df1$start, end = df1$end), seqnames=df1$chr)
df2.gr = GRanges(IRanges(start=df2$posn, end = df2$posn), seqnames = df2$chr)

df1.git  = GNCList(df1.gr)
t1 = Sys.time ()
overlap_git = findOverlaps(df2.gr, df1.git)
t2 = Sys.time()
sub1 = difftime(t2, t1, tz, units = c("auto"))
print (sub1)
print (overlap_git)

t1 = Sys.time ()
overlap_intersect = intersect(df1.gr, df2.gr)
t2 = Sys.time ()
print(overlap_intersect)
sub2 = difftime(t2, t1, tz, units = c("auto"))
print(sub2)

One should be informed about differences between IRanges and GRanges (for example chromosome context); for example, findOverlaps with IRanges returns intersections from both chr1 and chr12, but with GRanges it doesn't
Interval Tree based search is way faster (almost 10x) than simple intersect function


Merging files gives memory allocation error I have two files like

> head(a)
  chr position
1   1   136962
2   1   562020
3   1   672948
> class(a)
[1] "data.frame"
>

And

> dim(b)
[1] 1855235      10
> head(b[,1:3])
            SNP chr position
1:  rs62635286          1    13116
2:  rs75454623          1    14930
3:    rs806731          1    30923
4: rs200943160          1    49298
5: rs116400033          1    51479
6: rs141149254          1    54490
> class(b)
[1] "data.table" "data.frame"
>

I want to find corresponding SNP in file b matched to position in file a so I did like this that gives memory allocation error.

> unique(merge(a,b,x.by=a$position,y.by=b$position))
Error: cannot allocate vector of size 137.2 Gb

I also tried rstudio cloud but uploading the bigger file never completes

Can you help me with task please?

The expected output would be something like

rs62635286  1   13116
rs376723915 1   91515
rs147061536 1   92858
rs62642131  1   135982
rs371474651 1   158006
rs201293782 1   665401
rs138476838 1   668374
rs1401137   1   2063094
rs3128293   1   2065339

  EDITED

> setDT(a)
> setkey(a, chr, position)
> setkey(b, chr, position)
Error in setkeyv(x, cols, verbose = verbose, physical = physical) :
  some columns are not in the data.table: chr,position
> setkey(b, Chromosome, Position)
> a[b, nomatch=0]
Error in bmerge(i, x, leftcols, rightcols, roll, rollends, nomatch, mult,  :
  Incompatible join types: x.chr (factor) and i.Chromosome (integer). Factor columns must join to factor or character columns.
>

You could try something like this :

First you should work with the same class of object : data.table for fast merging operation

# Turning a to a data.table
setDT(a)
# b is already a data.table

Then you might have to set the keys of both table :

# I added chr in case you need to match the same position on a different chromosome.
setkey(a, chr, position)
setkey(b, chr, position)

Finally perform the joining operation :

a[b, nomatch=0]

Because you didn't provide the expected output, it might need some tuning !


Calculate the occurrence of motif in sequences (including overlaps) using Python I need to calculate the occurrences of a motif (including overlaps) in sequences (motif is passed in the first line of standard input and FASTA sequences in subsequent lines). The sequence name starts with >, and after whitespace is just a comment about the sequence that needs to be neglected. The input of program is like:

 AT
 >seq1 Comment......
 AGGTATA
 TGGCGCC
 >seq2 Comment.....
 GGCCGGCGC
 GCCGGCGGG

The output should be:

seq1: 2
seq2: 0

I decided to save the first line as a motif, strip the comment from sequence name, join lines of sequence in one line and save sequence names (keys) and sequences (values) in a dictionary. I also wrote a function for motif_count and want to call it on dictionary values and then save it in a separate dictionary for final output. However, it does not work the way I want it as the method motif_count does not produce any result when called on my d values.

Here is my code:

 #!/usr/bin/env python3

 import sys

 sequence = sys.stdin.readlines()
 d = {}
 temp_genename = None
 temp_sequence = None
 motif = sequence[0]

def motif_count(motif, temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
    if start <0:
        return count
    count += 1

if sequence[1][0] != '>':

    print("ERROR")

exit(1)

for line in sequence[1:]:
    if line[0] == '>':
       temp_genename = line.split(' ')[0].strip()
       tmp_genename = temp_genename.replace('>', '')
       temp_sequence = ""

else:

      temp_sequence += line.strip()

d[tmp_genename] = temp_sequence
for key in d:
    motif_count(motif, key)

This is a good task for Biopython:

from Bio import SeqIO

with open('input_file') as f:
    motif = next(f).rstrip().lower()
    for r in SeqIO.parse(f, 'fasta'):
        print(r.id, ':', r.seq.lower().count_overlap(motif))

Note: in future I would recommend don't add the motif to the start of your FASTA file, keep it separate so as not to alter the standard file format.


How to extract the protein fasta file from a genbank file? Let's say I have a genbank file, e.g. the one downloaded from here, that contains entries as

 gene            190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
 CDS             190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /function="leader; Amino acid biosynthesis: Threonine"
                 /GO_process="GO:0009088 - threonine biosynthetic process"
                 /codon_start=1
                 /transl_table=11
                 /product="thr operon leader peptide"
                 /protein_id="NP_414542.1"
                 /db_xref="ASAP:ABE-0000006"
                 /db_xref="UniProtKB/Swiss-Prot:P0AD86"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
                 /translation="MKRISTTITTTITITTGNGAG"

and I want to extract the potein fasta from it, how would I do this?

I found a code sample here. When I then run:

from Bio import SeqIO

gbk_filename = 'NC_000913.3.gb'
faa_filename = 'NC_000913.3_warwick.fasta'
input_handle  = open(gbk_filename, "r")
output_handle = open(faa_filename, "w")

for seq_record in SeqIO.parse(input_handle, "genbank") :
    print "Dealing with GenBank record %s" % seq_record.id
    for seq_feature in seq_record.features :
        if seq_feature.type=="CDS" :
            assert len(seq_feature.qualifiers['translation'])==1
            output_handle.write(">%s from %s\n%s\n" % (
                   seq_feature.qualifiers['locus_tag'][0],
                   seq_record.name,
                   seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()

it results in

     31     for seq_feature in seq_record.features :
     32         if seq_feature.type=="CDS" :
---> 33             assert len(seq_feature.qualifiers['translation'])==1
     34             output_handle.write(">%s from %s\n%s\n" % (
     35                    seq_feature.qualifiers['locus_tag'][0],

KeyError: 'translation'

How could this be fixed?

One can get it to work by using SeqIO.InsdcIO.GenBankCdsFeatureIterator:

from Bio import SeqIO

file_name = 'NC_000913.3.gb'

# stores all the CDS entries
all_entries = []

with open(file_name, 'r') as GBFile:

    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)

# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

This will then give the desired output (just the first two entries):

>b0001
MKRISTTITTTITITTGNGAG
>b0002
MRVLKFGGTSVANAERFLRVADILESNARQGQVATVLSAPAKITNHLVAMIEKTISGQDA
LPNISDAERIFAELLTGLAAAQPGFPLAQLKTFVDQEFAQIKHVLHGISLLGQCPDSINA
ALICRGEKMSIAIMAGVLEARGHNVTVIDPVEKLLAVGHYLESTVDIAESTRRIAASRIP
ADHMVLMAGFTAGNEKGELVVLGRNGSDYSAAVLAACLRADCCEIWTDVDGVYTCDPRQV
PDARLLKSMSYQEAMELSYFGAKVLHPRTITPIAQFQIPCLIKNTGNPQAPGTLIGASRD
EDELPVKGISNLNNMAMFSVSGPGMKGMVGMAARVFAAMSRARISVVLITQSSSEYSISF
CVPQSDCVRAERAMQEEFYLELKEGLLEPLAVTERLAIISVVGDGMRTLRGISAKFFAAL
ARANINIVAIAQGSSERSISVVVNNDDATTGVRVTHQMLFNTDQVIEVFVIGVGGVGGAL
LEQLKRQQSWLKNKHIDLRVCGVANSKALLTNVHGLNLENWQEELAQAKEPFNLGRLIRL
VKEYHLLNPVIVDCTSSQAVADQYADFLREGFHVVTPNKKANTSSMDYYHQLRYAAEKSR
RKFLYDTNVGAGLPVIENLQNLLNAGDELMKFSGILSGSLSYIFGKLDEGMSFSEATTLA
REMGYTEPDPRDDLSGMDVARKLLILARETGRELELADIEIEPVLPAEFNAEGDVAAFMA
NLSQLDDLFAARVAKARDEGKVLRYVGNIDEDGVCRVKIAEVDGNDPLFKVKNGENALAF
YSHYYQPLPLVLRGYGAGNDVTAAGVFADLLRTLSWKLGV


predict the IC50 after using Schrodinger suite Glide tool for docking How to predict the half-maximal inhibitory concentration (IC50) of my ligands.
I'm currently using Schrodinger suite software and I don't know how to predict the IC50?

I don't think it is possible just from the docking scores. Docking scores are designed to give you an estimate of the binding affinity (entalphy), but they do not take entropic contributions into account very well. The best you can do are molecular dynamics simulations to calculate the free energy of binding. Which can be computationally quite costly e.g. with GROMACS or CHARMM.
You can estimate the IC50 with molecular dynamics simulations as well, but that is not straightforward. E.g. with Umbrella simulations. An example for a potassium channel and a very simple ligand can be found here.
Another approach would be to build a predictive model. If you have a library of compounds at hand you could use the chemical structures and e.g. the docking scores to estimate the IC50 based on statistics or machine learning. But biology is very messy in general, though if you get a representative dataset with very similar compounds to the one that you have in your dataset then you might be successful. You have to be sure that the binding site and pose are correct though.


ATAC seq peak annotation : Which peak to choose for downstream analysis After doing peak call I'm annotating the peaks using chipseeker tool, which I want to take further downstream analysis. Now unlike RNA seq where i have single gene and their respective expression. But here the peaks for a single gene multiple peaks. As an example have filtered promoter peaks.
seqnames  start    end width strand     V4 annotation geneChr geneStart geneEnd geneLength geneStrand geneId
1     chr1 826797 828101  1305      *  Peak1   Promoter       1    826832  852225      25394          1 643837
2     chr1 869647 870324   678      *  Peak4   Promoter       1    868240  870201       1962          2 284593
3     chr1 876462 877140   679      *  Peak5   Promoter       1    868071  876903       8833          2 284593
4     chr1 877219 878105   887      *  Peak6   Promoter       1    874529  877234       2706          2 284593
5     chr1 923459 924318   860      * Peak14   Promoter       1    923928  939291      15364          1 148398
6     chr1 924379 926062  1684      * Peak15   Promoter       1    925150  935793      10644          1 148398
       transcriptId distanceToTSS         ENSEMBL    SYMBOL                                    GENENAME
1 ENST00000623808.3             0 ENSG00000228794 LINC01128 long intergenic non-protein coding RNA 1128
2 ENST00000432963.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
3 ENST00000446136.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
4 ENST00000635557.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
5 ENST00000420190.6             0 ENSG00000187634    SAMD11    sterile alpha motif domain containing 11
6 ENST00000437963.5             0 ENSG00000187634    SAMD11    sterile alpha motif domain containing 11

Now if we see first 6 rows we have two genes with three different peaks. When i load and see in the igv browser I do see the there are subtle differences in the accessibility signals.
How to address this issue ,in other words which peak to consider or take all the peak and consider their accessibility as one.
Any suggestion or help would be really appreciated.

You can see from your table that the different peaks correspond to different promoters for different transcripts, for example, FAM41C has ENST00000432963.1, ENST00000446136.1, ENST00000635557.1. And the fact that they are all within the TSS (see distanceToTSS) suggest these peaks are actually falling onto actually TSS.
The underlying issue is are these TSS or do will we mis classify the peaks that are inside genes or other non-coding elements as TSS peaks. Nothing to do with TF, histones etc.
I cannot see the number of reads in those peaks, but if one of them is really strong, and the rest almost non-existent it would make sense.
How can we improve this, I think it depends on the aim of your analysis. Since you are looking at TSS, I guess it makes sense to be able to link it to something that is actually expressed. You can do:

To use an updated annotation (i see that some of the transcripts come from old hg37 annotation)

you have prior knowledge which transcripts might be expressed from RNA seq data, I would use that to reduce the number of transcripts


How to retrieve logical expressions (KO based) for reactions from KEGG? The completeness of a module can easily be checked by looking at the Definition entry associated with the module. For example, in module M00010, it is given as:

Definition  K01647 (K01681,K01682) (K00031,K00030)

which can be translated to:

K01647 AND (K01681 OR K01682) AND (K00031 OR K00030)

If this expression evaluates to TRUE, the module is considered complete.

I am wondering whether analogous information exists for a single reaction. For example, when looking at reaction R00352, one finds the following information about its Orthology:

But I can't easily find the logical relation that correspond to the KOs for this reaction. It could be:

K01648 AND K15230 AND K15231

or:

K01648 OR K15230 OR K15231

or:

K01648 OR (K15230 AND K15231)

and so on.

In the above example, I know that the correct expression is:

K01648 OR (K15230 AND K15231)

i.e. one either needs either K01648 or both of the other two subunits.

Can this information be retrieved from KEGG for each reaction and if so, how?

Some information seems to be available on the module level utilising the API.

Reaction R00352 points to the module M00173. While some information is in the Definition, it's not trivial to link it to the particular reaction. This linkage is however easier accessible via API:


Counting letters in phylip alignment columns with Biopython I have been using python 3.6 and biopython 1.72 to work with protein data files. I am using a protein sequence file (phylip format), for example:

    14    678
Zebrafish  LSSCGVVSGD LISVILPASS LEETQTSSAA AHQTHTDQQA GGSHVSSSSS
Fugu       LASCGIVSGD LIRVALIRAA LMAADAPDRA DDGGGHSEQA VAMATNQVS-
Chicken    LASCGIVSGD LIRVALIRAA LMAADAPDRA DDGGGHSEQA VAMATNQVS-
Zebra      LASYGIVPGD LICLLLEETD AKSSLPPPSS SSPPSLQNGR ETSTLTPKNQ

I want to manipulate the file in a way so that;

I may get the number of columns which have all the same letters (amino acids)
The number of columns which have at least 2 different amino acids (with no dashes).
The number of columns which have any single dash.

I have iterated over the columns of the file with:

alignment = AlignIO.read(open("example.phy"), "phylip")
def column(alignment):
    n = len(alignment[0].seq)

    for i in range(n):
        c = []

        for record in alignment:
            c.append(record.seq[i])

        yield c

This returns the columns of the file, but I don't understand how I can manipulate them further.

Can I get some helpful suggestions?

using something like for col in column(alignment):

I may get the number of columns which have all the same letters (amino acids)

len(set(col)) == 1

The number of columns which have at least 2 different amino acids (with no dashes):

len(set(col)) > 1 and '-' not in col

The number of columns which have any single dash:

'-' in col

or if you really want columns with a single dash:

    col.count('-') == 1


Salmon output: transcripts quantified with zero reads support I quantified some samples using Salmon. According to the documentation of output format, the last column of the 'Quantification File' represents the number of reads that supported the given transcript, namely:

  NumReads — This is salmon’s estimate of the number of reads mapping to each transcript that was quantified. It is an “estimate” insofar as it is the expected number of reads that have originated from each transcript given the structure of the uniquely mapping and multi-mapping reads and the relative abundance estimates for each transcript.

In some case, I had in the output some transcripts that had NumReads=0.000000. Does this mean that I should discard them from the list of transcripts obtained in that sample? How comes that those transcripts end up in the output file in first place, since they don't have any reads support and the quantification was done using the reads themselves?

All transcripts should always be returned in the output, whether there's evidence supporting their expression or not. This is for the sake of simple convenience, since typically one would be merging these values across samples into a matrix. Missing values will then tend to become either NA or cause an error, when what would be appropriate would be to have values of 0. This also facilitates calculations such as "what's the median number of reads per transcript?" or "what percent of transcripts have an expected reads count of at least X".


co-localisation error code using co-loc Can anyone offer advice on what the below error code relates to. I am using co-loc R package with 2 GWAS datasets. Using coloc.abf function with pvalues, MAF, sample size and case/control cohort. Thank you

  results <- coloc.abf(  dataset1=list(snp=as.character(1_subset[,2]),pvalues=1_subset[,11],N=1_data_subset$sample_size, type="cc",s=0.327,MAF=1_data_subset[,12]),   dataset2=list(snp=as.character(2_subset[,2]),pvalues=2_subset[,8],N=2_subset$sample_size,s=0.4319,MAF=1_subset[,12], type="cc"),p12=1e-06)

error code: Error in data.frame(V, z, r, lABF) :
  arguments imply differing number of rows: 1994, 1992
In addition: Warning messages:
1: In 2 * N * f :
  longer object length is not a multiple of shorter object length
2: In 2 * N * f * (1 - f) :
  longer object length is not a multiple of shorter object length
3: In r * z^2 :
  longer object length is not a multiple of shorter object length

In your dataset2 = list(), the MAF argument should be MAF=2_subset[,12] instead of MAF=1_subset[,12]. Since your two datasets have different numbers of lines, this is what causes the error.

I am not sure whether this is a typo in your code, or it is because you do not have MAF information in the second dataset. What I like to do before COLOC is to filter two datasets so they have the same number of lines, and then feed those filtered datasets to the coloc function.

(Also just out of curiosity, is 1_subset a data frame for your GWAS results? I think it is not recommended to have variable names starting with a number.)


Error in adding data in columns of R data frame This is my file

df1 = read.csv("HSC_LSC_BLAST_karyotyope.txt",header = TRUE,sep = ",",row.names = 1)

metadata <- data.frame(row.names = colnames(df1))
metadata$Group <- rep(NA, ncol(df1))
metadata$Group[seq(1,4,1)] <- 'HSC'
metadata$Group[seq(5,15,1)] <- 'Blast'
metadata$Group[seq(16,23,1)] <- 'LSC'
#metadata$Group[seq(13,16,1)] <- 'Mono'
metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy'
metadata$Cytogenetics[seq(5,21,1)] <- 'NK'
metadata$Cytogenetics[seq(22,23,1)] <- 'Abnormal'

>s metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy' Error in `$<-.data.frame`(`*tmp*`, Cytogenetics, value = c("Healthy",  :
> replacement has 4 rows, data has 23

Im trying to add other metadata information based on clinical information apart from giving sample labeling , im certainly doing something wrong but im not sure what it is. As my idea is the another column named Cytogenetics which would have respective information status but im getting that error.

Add the Cytogenetics columns and then assign things to a subset of it.

metadata <- data.frame(row.names = colnames(df1))
metadata$Group <- NA     # you don't need to write a vector; just NA does the job
metadata$Group[seq(1,4,1)] <- 'HSC'
metadata$Group[seq(5,15,1)] <- 'Blast'
metadata$Group[seq(16,23,1)] <- 'LSC'
metadata$Cytogenetics <- NA    # initiate the column with NA values
metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy'
metadata$Cytogenetics[seq(5,21,1)] <- 'NK'
metadata$Cytogenetics[seq(22,23,1)] <- 'Abnormal'


Rearrange your text file into a sorted BED file:

$ tail -n+2 segmentation.txt | awk -vFS="\t" -vOFS="\t" '{ print $2, $3, $4, $1 }' | sort-bed - > segmentation.bed

Use BEDOPS bedmap with --echo-map-id-uniq and --count operators on merged regions:

$ bedmap --delim '\t' --echo --count --echo-map-id-uniq --fraction-map 1 <(bedops --merge segmentation.bed) segmentation.bed > answer.bed

The file answer.bed will contain merged intervals from segmentation.bed in the first three columns, the number of intervals from segmentation.bed which are contained entirely within the merged interval in the fourth column, and a unique listing of IDs of those contained intervals in the fifth column.


DeSeqDataSet experimental design: column with integer values I am creating DeSeqDataSet and I am unsure whether I need to include the number of mice into the design formula. Basically, the samples look like that:

I have 4 different groups: GF, B6, ET and Wld. It is obvious that these should go into the design formula, but I am not sure about the mice number. Basically for each group, e.g. GF, there are 4 experiments with 1, 2, 3 and 4 mice each. When I am trying to use both Groups and Mouse columns as a design formula, R is telling that I need to convert Mouse to factor:

I am a bit confused about two things here: first, as I told already, should I include this Mouse column at all in the design formula of DeSeqDataSet and second, if I do need it, how to convert it to factor? Any suggestions would be greatly appreciated.

  Update

What concerns me also is that if I specify design formula as ~ Groups + Mouse, the DeSeq() function starts comparing the Mouse but not the Groups which is the most interesting for me.

Any IDs (e.g. replicate number) don't need to be included in the design formula if it's not a repeated factor in multiple samples.

In this case, you have said that mouse number is related to the number of mice used in the pool (and is a variable factor that is consistent across multiple samples), in which case it should be included in the design matrix. This would also be the case if the same mouse were being used for multiple RNASeq experiments. Converting an integer value to a factor is done by the factor method:

sampleTable$MiceUsed <- factor(sampleTable$MiceUsed);

Other examples of things that might be useful in the design formula would be sex, age, and sequencing run batch (assuming these are not consistent across samples).

With a design formula of something like ~Groups + Mouse, it will also be comparing Groups, it's just hidden in the results. The DESeq2 vignette is really great at explaining what's going on, but as a quick answer you can pull out the Groups results with the contrast function. This will pull out group differences that are independent of mouse number:

results(dds, contrast=c("Groups", "GF", "Wld"))

[first value is the variable name, second value is the numerator condition, third is the denominator condition]


GATK 4.1.4.0 Mutect stats output I have adapted my pipeline to the new filtering strategies implemented in gatk 4.1.4.0 and while looking at its output I noticed that the stat file generated by mutect2 has a negative number:

callable -1.538687311E9

I did not find the documentation about this file content but in the gatk forum it is mentioned that is should represent the number of bases above a certain coverage (10 by default).
This happened for a WGS sample, while for WES I obtain a more sensible value (well, at least positive and integer, even if maybe a tad high):

callable        8.6884303E7

Has this happened to anyone else?
I compared the resulting calls obtained for the WGS with a previous version of mutect2 and they are comparable-ish (98% of the new "ones" are in common with the "old" ones), but I am still worried that something fishy is happening, expecially considering that the filtering step of the pipeline uses this stats file.

This was answered by the OP -- it was a bug that has been fixed in a subsequent release, 4.1.5.0. See the release notes here: https://github.com/broadinstitute/gatk/releases/tag/4.1.5.0

  Made callable sites a Long to avoid integer overflow (#6303)

Issue ticket: https://github.com/broadinstitute/gatk/issues/6302


KEGG FTP vs KEGG API I was reading the KEGG plea and I found that it doesn't forbid using the KEGG API. Then, what is in the FTP server license for personal use/academic use that it is not covered by the API?

Or I could download all the database via the API?

PS: I couldn't find how to download which genes are in which pathways but it might be I didn't expend enough time with the queries.

I found the way to download the genes and pathways thanks to this answer using KEGGREST.

It seems that the API allows to fully download the database because it allows to download info, list, find, get, conv, link, ddi.

In fact in the message of 2011-05-21 it mentions that the API includes more options to download the KGML files that previously where only possible to get via the FTP site.


Looking for a taxonomic dataset that links tree genera (or species) to families I am looking for a taxonomic dataset that would allow me to links tree genera (or species) to families. In particular, I am interested in the tropical species found in the Amazon.

I found a long but not exhaustive list on Wikipedia:
https://en.wikipedia.org/wiki/List_of_trees_and_shrubs_by_taxonomic_family

After processing it, I realised it contains information about 404 genera, while the dataset I am dealing with has around 600. (in the end, I only got 23% of my data covered, very sad)

Does anyone know if such a crosswalk file exists, in some scientific article or somewhere in the Internet?

Moreover, please tell me if there are other StackExchange websites which could be more appropriate, I wasn't sure of where to post this.

The NCBI taxonomy is a quite extensive resource: https://www.ncbi.nlm.nih.gov/taxonomy

I would thus expect there is a good chance that it covers all of the taxa that you are interested in.  It officially claims that it is not an authoritative source, but it's still pretty darned good.


Create a NHX phylogeny format with r or python I need help to add information in my phylogenetic tree on each branch.
I have two pieces of information at my disposal:
1 ) A phylogenetic tree in newick format :
**(mammal:0.14,turtle:0.02,(rayfinfish:0.25,(frog:0.01,salamander:0.01)50:0.12)95:0.09);**

2 ) A matrix with the edge number and information in the Information column:
Edge    Inf
7        0
6        0
5        2
4        1
3        2
2        2
1        2

Does someone have an idea in order to add these informations by creating a NHX format in python (witg ete3 for instance) or with a R package ?
Here is the tree :

Thank you for your help

The treedata book has a section on tree annotations that should work for you.


How can I work out gram staining information from the species name? I have a very long list of prokaryote organism names. I would like to know their gram staining (gram positive or gram negative) information.

Please let me know if there is any database that keeps this information.

I am also trying to find a program that can find gram staining information via organism name only as my organism list is long.

Edit: for each species the NCBI taxonomy lineage is also available

The GOLD database can be used to retrieve a list of species names according to their Gram staining: go to the "Advanced search" page and select the either "gram-" or "gram+" under "Organism Fields" -> "Gram Stain". Then on the results page click on the number under "Organisms". On the results page you can select to have 1000 results per page, meaning that you would need to do some manual work.

I should add that as far as I know this information is based on what people submit when they register their genome projects, so I wouldn't expect this list to be curated or 100% accurate.

If you have access to the genome of the species you are interested in (you could download a representative strain for each species using ncbi-genome-download), you could use traitar, which can predict the gram staining, among other interesting broad phenotypes.


Dataset: Locations of regulatory sequences in the human genome? I am looking for the positions of annotated regulatory sequences (promoters, enhancers and suppressors) in the human genome.

I looked at Ensembl regulatory Build and PAZAR but I am not used to look for datasets and I failed to find what I was looking for.

Can you help to find annotated regulatory sequences (start and end positions) in the human genome?

This can be done quite easily using Ensebl's BioMart.

Choose the Ensembl Regulation database:

Select the "Human Regulatory Features" dataset:

That's basically it right there, just click on "Results":

Export to file and click "Go":

This will download a file called mart_export.txt which looks like this (I chose TSV for tab separated values):

$ head mart_export.txt
Chromosome/scaffold name    Start (bp)  End (bp)    Feature type
18  76429380    76430144    Open chromatin
8   66405962    66406502    Open chromatin
4   61184401    61184600    CTCF Binding Site
X   40733600    40737000    Promoter
5   97407001    97407200    CTCF Binding Site
X   73946201    73946600    Promoter Flanking Region
15  19948201    19949200    CTCF Binding Site
5   11302601    11303143    Open chromatin
2   208407801   208408000   CTCF Binding Site

As you can see the format is pretty straightforward. The first column has the chromosome or scaffold name (it should always be a chromosome for human), the second and third are the star and end positions and the last field is the type of region.

You can also choose to limit by region type by choosing a filter (click the "Filters" link) in BioMart before downloading or by simply parsing the file once you've downloaded it.


How do I interpret the function O(f^-2 log n) in Shekhar et al. (2017)? The results of Shekhar et al. (2017) "indicate that ASTRAL requires

$\mathcal{O}(f^{-2} \log{n})$

gene trees to reconstruct the species tree correctly with high probability where n is the number of species and f is the length of the shortest branch in the species tree.

How do I interpret this function, specifically the big O notation, when applying it to my data to estimate whether I have enough loci?

Big O notation means "of the order of". Thus $\mathcal{O}(f^{-2} \log{n})$ means that the number of trees required goes with $\log{n}$ for number of species (i.e. if you square the number of genes, you should double the number of trees) and with the inverse square of the shortest branch length (i.e. if you double the shortest branch length, you quarter the number of trees needed).

Big O notation generally allows to to compare two input sizes, but cannot really tell you anything about a single point. I.e. it is hard (impossible?) to tell from this function whether you have sufficient loci.


Duplicated genes and genome fragmentation This is a follow up question to this post: How can I find out if my gene of interest is duplicated?

Briefly, to identify homologous genes, I performed tBLASTn analyses against my genome of interest using various sequences of a closely related species as queries for the BLAST analyses. Using this method, I found two (and in some cases up to five copies) of my different genes of interest in the target genome. These genes are involved in different organelle functions (i.e., they are not surface antigens) and I want to find out whether these sequences represent duplication events or not. These analyses are being performed on protozoans parasites that belong to the same genus. The query genome is properly annotated, but the target genome is not (it is a whole genome shotgun project that has been assembled at the scaffold level).

Following the advice of @Michael G., I re-analyzed my sequences and blasted them against the entirety of Genbank. If there are two sequences for the same protein (for e.g. protein 3-1 and protein 3-2) and I compare them among themselves they are very similar (over 70% identity at the protein level). When I blast each sequence against the whole protein RefSeq database, the first hit is my query and the rest of the hits correspond to the same protein in other organisms. I believe that these sequences represent duplication events and it is possible that they are paralogues, but more evidence is needed to classify them as such.

However, when there are more than two sequences for the same protein, I found something different. For e.g., say I have three sequences for the same protein: sequence 5-1, sequence 5-2 and sequence 5-3. Sequences 5-2 and 5-3 are very similar (78% identity at the protein level) but when I compare sequences 5-1 and 5-2 I find that they are overlap in the middle, i.e., they are part of the same sequence. Thus, it appears that I don't really have three sequences for this protein but two. Blastp analyses against the whole RefSeq database give the expected result (first hit is my query sequence, the rest of the hits are same sequence in other organisms). So, once again I have what appears to be a duplication event, but in this case one of the sequences appears to be split in two -- could this be because the target genome is not completely assembled? Why could this be happening?

Thank you very much for your help.

-Leah

Edit: The new results are posted in this thread: Duplicated genes and genome fragmentation - Part 2

The tBlastn looks good and looks like some sort of parology, although the hits you are getting need to be checked to make sure they are not psuedo-genes. The shot-gun genome however appears poorly annotated because it had misannotated 5-1 and 5-2. What I don't understand is why you haven't performed a blast back over the original highly annotated genome? You must have done this?

There are a few avenues to go down here, what I'd like to know which group of protists are you working with? I don't need to know the species, but I do need to know the family please.


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

MUMmer includes a repeat-match tool which can do exactly that:
$ repeat-match
USAGE:  repeat-match  [options]  <genome-file>

Find all maximal exact matches in <genome-file>

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

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


Highly variable genes analysis by DESeq2 In the following article. There are two ways of getting highly variable genes:

resSig <- subset(res, padj < 0.1)
head(resSig[ order(resSig$log2FoldChange, decreasing=TRUE), ])

Getting genes that are statistically significant, then ordering them by the log2FoldChange which is

  The column log2FoldChange is the effect size estimate. It tells us how much the gene’s expression seems to have changed due to treatment with dexamethasone in comparison to untreated samples.

And the second one:

topVarGenes <- head(order(rowVars(assay(rld)),decreasing=TRUE),20)

I am struggling to understand what are the differences between the two ways of getting the most variable genes. What is the difference?

The first method (resSig$log2FoldChange) is for getting the most differentially expressed genes. These are the genes with the biggest differences between specific groups you pre-define.

The second method (rowVars(assay(rld))) is for getting the highly variable genes. These are the genes most variable across all samples regardless of which samples they are.


Finding common and unique data set by comparing two files based on their column and to split the columns multiple strings to print in output I have  very large sizes tab-delimited .vcf files and want to match these two / or 3 files based on their position and print to a new .csv file

File structures:

File_1: tab-delimited file (.vcf) and its as column names as follows

(line number 3439) #CHROM     POS     ID       REF      ALT      QUAL    FILTER
 INFO

File 2: same as file_1 column names

(line number 3407) #CHROM     POS     ID      REF     ALT     QUAL    FILTER
 INFO  FORMAT  SAMPLE_1  SAMPLE_2 .....

In file_2, column 7 (INFO) contain many substrings like
AC=46,2;AF=0.958,0.042;AN=48;DP=269;ExcessHet=3.0103;FS=0.000;InbreedingCoeff=0.5411;MLEAC=92,4;MLEAF=1.00,0.083;QD=25.36;SOR=2.488  from these strings , have to print  only information of AC=, AF=, AN=, DP=

Desired output files to generate:

position matched in two of files: common_position_matched.csv

if I have more than 3 files also, the output should be one file and important thing, if a position column (1-POS) only matched only in 2 files and 1st file line should be NA NA NA

file_1  CHROM  POS    REF  ALT     file_2 #CHROM  POS      REF  ALT  INFO

1       22    10511521  C     T         1    22  10511521   C    T   AC=46,2 AF=0.958,0.042 AN=48 DP=269
2       22   10510544   G     A         2    22  10510544   G    A   AC=49,2 AF=0.958,0.042 AN=89 DP=536
3       22   10515068  AGAT,T AGAT,AT   3    22  10515068   AAA  AAAGG,A,GAA AC=100  AF=0.958,0.042 AN=62 DP=123
4       22   10515118  A G,   TAA       4    22  10515118   AG,  TAA AC=32   AF=0.958,0.042 AN=45 DP=500
5       22   10515118  AAAG   A         5    22  10515118   AATG A   AC=50   AF=0.958,0.042 AN=49 DP=129

note: while doing matching, not removing the duplicates, because in the same position there may be an addition or sometimes it may be deletion.

unique position of each file, in tab-delimited

output: File1_unique.csv and File2_unique.csv etc.

so far was able to read the file and match them according to position and print the output, but was not able to write efficient code

import pandas as pd
df1 = pd.read_csv("File1_3.vcf",sep='\t',usecols = ['POS']) ## Reading file1
df2 = pd.read_csv("file2_3.vcf", sep="\t", usecols = ['POS']) ## Reading file2
df3 = pd.concat([df1,df2], sort=True) ## Combining both the dataframes
df4 = df3.drop_duplicates(keep=False) ## Dropping the duplicates (intersect)
df4.to_csv("c3-UniquePosition_of_bothData.csv", sep="\t", index=False, header=True) ## Writing the unique to both
df1_Uni_file1_c3Posi = pd.merge(df4, df1, on='POS', how='inner') ## Identifying the unique position of File1
df2_Uni_File2_c3Posi = pd.merge(df4, df2, on='POS', how='inner') ## Identifying the unique position of File2
df_File1_File2_common_c3Posi = pd.merge(df1, df2, on='POS', how='inner') # Identifying the common chr-position of File1 and File2```

Program 2: (giving original file without editting)
import pandas as pd
df1= pd.read_csv("File1_22.vcf.gz", sep="\t", skiprows=3438, usecols = [0,1,2,3,4])
df2 = pd.read_csv("File2_22.vcf.gz", sep="\t", skiprows=3406, usecols = [0,1,3,4,7])
#writing the output files
#df1.to_csv("File1_c22.csv")
#df2.to_csv("File2_c22.csv")
#mergeing
df3 = pd.merge(df1, df2, on='POS', how='inner', sort=True)
df3.to_csv("common_position.csv", sep=",", index=False, header=True)
#df3 = pd.concat([df1,df2], axis=1).to_csv('check1.csv') # this command join multiple output to single output

Could any one give efficient python pandas script, to do this

Thanks All

One way to do this is to use a dictionary. You can read the first file line by line, and for each line, use "chr:position" as a key, and the corresponding line as the value. If you want to handle duplicates, you can make the value a list and store multiple lines with duplicated "chr:position" as elements in that list.

After constructing this dictionary, read the second file line by line. For each line, set "chr:position" as the so-called name of that line. Then using this "chr:position", you can easily find the corresponding line in the dictionary. Append the line from the second file to that line from the first file.

After finishing reading the second file, you are ready to write the new file to output. To do so, iterate over the keys in the dictionary, and write the value to each key to a new file. Remember that each value in that dictionary is now a new line.

I choose to use a dictionary because they are basically a hash table, meaning that using the key you can instantly read that value in the dictionary, instead of iterating every element in a list to find that element. This is very handy when you have a large dictionary.


What kind of analysis can be done with differential expression of transcription factors? I have two different stem cell types and their respective gene expressions from RNAseq.

I noticed that there is differential expression in some of the genes coding for histones as well as transcription factors?

What kind of interesting downstream analysis might be available?

Is it feasible to survey ENCODE for downstream TF binding sites?
Can anything interesting be investigated with regard to the histone gene data?

Thanks in advance.

Expanding my comment.

TF usually bind in the upper region of the genes it regulates, so you could see which genes are in the region below the binding site. Perhaps there are some databases for TF.
I found a website base on the ENCODE project with the genes regulated by each TF.You can use them for pathway or gene set enrichment analysis, but I am not aware of any specific analysis for genes regulated by a single TF.

The histone data is harder if it was methylation I would say to look if different states of the histones are in different spatial state, but with gene expression is harder to conclude anything. Nevertheless you can check methylation studies to see if gene expression is associated with methylation changes.


What is the "primer masterfile" format? I am working with an amplicon kit from Swift Bioscience and they provide the primer locations in what they call a "masterfile". The format is like this:

$ head -n2 Acel-Amplicon_56g_masterfile.txt
1       43814951        43815064        MPL_S505N-A519T_F/2     43814932        43814951        MPL_S505N-A519T_F/2_primer      43815064  43815081 MPL_S505N-A519T_R/1_primer      GCCTGGATCTCCTTGGTGA     CTTCGGCTCCACCTGGT
1       115256480       115256576       NRAS_G60-S65_new_R/1    115256457       115256480       NRAS_G60-S65_new_R/1_primer     115256576 115256602        NRAS_G60-S65_new_F/2_primer     TGATGGCAAATACACAGAGGAAG CAGGATTCTTACAGAAAACAAGTGGT

This is the format which is understood by Swift's primer clipping tool, but in order to use other tools, I would like to convert it to BEDPE so I can then use standard tools to clip the primers from my bam files. However, I can't do that since I am not sure what the fields in this file are.

This is what I think thy are:

Chromosome
Target region start
Target region end
5' Primer name
5' Primer start
5' Primer end
3' Primer name
3' Primer start
3' Primer end
5' Primer sequence
3' primer sequence

Is my interpretation correct? Also, is this a standard file format or something specific to the Swift tool? I couldn't find any documentation about a "masterfile" format online.

this is the swift standard primer file. You are a bit off in your interpretation of the columns (seems like you missed one as there are 12 columns, not 11 which threw everything off):

1) Chromosome

2) Target Start

3) Target Stop

4) Target Name (change which affected everything downstream in counting)

5) 5' Primer Start

6) 5' Primer Stop

7) 5' Primer Name

8) 3' Primer Start

9) 3' Primer Stop

10) 3' Primer Name

11) 5' Primer Sequence

12) 3' Primer Sequence

To convert to a BEDPE describing the primer locations, you only will need columns 1, 5, 6, 8, and 9 (and optionally 4). The last letter of the Region name for the denotes the strand orientation of the region. For example Region 1 has _F denoting forward or (+) strand for BEDPE usage, while region 2 has _R in name denoting reverse or (-) strand for BEDPE. Putting all this together, you can convert to BEDPE with:

$ awk -v OFS="\t" '{
            o=".";
            if($4~/_F\//){
                o="+"
            }
            else if($4~/_R\//){
                o="-"
            }
            print $1,$5,$6,$1,$8,$9,$4,".",o}' foo.masterfile
1   43814932    43814951    1   43815064    43815081    MPL_S505N-A519T_F/2 .   +
1   115256457   115256480   1   115256576   115256602   NRAS_G60-S65_new_R/1    .   -


How to export peptide-shaker data from commandline? I would like to export the data from Protein-Shaker proteomics analysis software and use it in downstream analysis. For example the overview data, numbers of validated peptides per protein, confidence, etc. Ideally I would like to generate this output from commandline without the gui and, if possible to plain txt or csv.

You have to start PeptideShaker calling java with the -cp argument an input project file and provide a comma/space-separated list of integers for the reports you want to export:

java -cp PeptideShaker-1.16.37.jar eu.isas.peptideshaker.cmd.ReportCLI -in input.cpsx -reports "0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10"`

Here is the list of reports you can export:

-reports                   Comma separated list of types of report to export.
                           0: Certificate of Analysis,
                           1: Default Hierarchical Report,
                           2: Default PSM Phosphorylation Report,
                           3: Default PSM Report,
                           4: Default PSM Report with non-validated matches,
                           5: Default Peptide Phosphorylation Report,
                           6: Default Peptide Report,
                           7: Default Peptide Report with non-validated matches,
                           8: Default Protein Phosphorylation Report,
                           9: Default Protein Report,
                           10: Default Protein Report with non-validated matches,
                           11: Extended PSM Report,
                           12-n: Your own custom reports.

I don't know whether you can also use the SearchGUI output zip to generate this files. The program accepts zip files as input, but I think that means PeptideShaker project that are exported to zip files, but not SearchGUI output zips. So, in the end you might need to run the GUI at least one time, unfortunately.


bcftools filtering all files in a directory Probably a silly oversight on my part, but I'm trying to filter all the vcfs in a directory with bcftools using a simple loop. My basic command is working fine:

bcftools filter -i 'QUAL > 1000' -o filter/file1out.vcf file1in.vcf

but when I try to loop it, it echoes through each filename, but generates no output files or error messages.

for f in *.vcf; do echo "filtering $f"; bcftools filter -i 'QUAL > 1000' -o filter/$f_out.vcf $f; done

Am I not handling the output direction properly? Any help much appreciated!

I assume your shell is looking for the variable $f_out, but can't find it because you did not define it. It probably generated a hidden file with the name .vcf in your output directory.

Instead, use the following: -o ${f}_out.vcf. This will "protect/limit" the bash variable to just $f rather than $f_out.

When encountering issues like this it's worth putting echo in your loop right before your command, in this case bcftools (adding quotes whenever necessary). That would show you the command as how it is going to be executed and would show you, in this case, that the output name is incorrect.


DiffBind, diffferentially binding site I have data for 3 histone marks (2 for silencing and 1 for activation) each mark has three replicates.
when I run the diffBind package I have three contrast:
Silence1-Activation
Silence1-Silence2
Silence2-Activation

When I run dba.plotBox for any of contrasts I will have the same distribution of reads.like below:

I even tried to start with 1000 peaks of each one that showed higher fold change but again got the same pattern for plotbox. How can I see which peaks belong to which histone marks?

I got the answer from Rory Stark. the samples should have a unique ID. (https://support.bioconductor.org/p/132502/#132514)


Pattern mining from a genomic sequence I need to find the following pattern from a genomic sequence

TTCGB-N15-TTCGB
B represents C, G, or T
N15 represents A, T, G or C and these four base combination should have the length of 15 bases.

To do the same, I tried grep command, but I could not specify B and N15 in the grep command line. Therefore, please help me to find the same.

For example, I have a nucleotide sequences as follows

>ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT

Expected outcome,

  TTCGTGCATGCATATGCGCATTCGC

Thanks in advance.

echo ">ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT" | grep -oP  'TTCG[CGT][ATGC]{1,15}TTCG[CGT]'

[CGT] looks for one occurrence of one of the three.
[ATGC]{1,15} looks for "up to 15 letter long" combinations of the four bases.

Moreover, -P makes your regex call Perl-like and is required to make this regex work. Different regex flavors behave (slightly) differently.


Challenging benchmarks for supervised learning on sparse scRNA-seq data One challenging aspect of modeling scRNA-seq data is data sparsity, that is, scRNA-seq measurements typically suffer from large fractions of observed zeros (i.e. dropouts), where a given gene in a given cell has no unique molecular identifiers or reads mapping to it [1].
What would be an appropriate dataset to benchmark supervised learning algorithms that operate on single-cell RNA-seq data? In particular, I am looking for challenging datasets (preferably with labels, i.e. cell type classification) where existing methods struggle due to the sparsity of the data. Additionally, does it ever happen that the dropouts are confounders of the labels? (i.e. for a given cell type A the dropout probability of a gene X is higher than for another cell type B).
Any pointers are welcome.

You could try the PBMC 3k data from the Satija lab:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
For well-annotated data, there's the Single-cell proteo-genomic reference map:
https://www.biorxiv.org/content/10.1101/2021.03.18.435922v1
Data for that are available here.
If these datasets are too large (i.e. too good), you can use the existing counts as weighting for a random removal of reads. Genes with high expression are likely to be present regardless of sequencing depth, hence the count-adjusted removal.
FWIW, I don't think that supervised learning is a good fit for dealing with data sparsity in cDNA read counts. There are plenty of existing attempts to solve the issues by adjusting the assumed models, including research by the authors of commonly-used bulk differential expression tools (e.g. here).


Why Illumina if PacBio provides longer and better reads? PacBio provides longer read length than Illumina's short-length reads. Longer reads offer better opportunity for genome assembly, structural variant calling. It is not worse than short reads for calling SNP/indels, quantifying transcripts. Sounds like PacBio can do whatever Illumina platform can offer.

Would that be any reason for going with Illumina if PacBio's long reads can do everything and more?

There are so many reasons why one might want to prefer Illumina over PacBio (also note that it's a false dichotomy, at least Oxford Nanopore is a competitive sequencing platform):

The first (IMHO and the most common reason) is still the cost of both sequencing and the instruments. Illumina can sequence a Gbp of data for \$7 - \$93. PacBio sequencing is according the same webpage \$115 per Gbp, however at our sequencing center it's ~$200. Though ONT might have already a cheaper solution. Edit, I just found a google sheat with prices that seems to be frequently updated, seems that the ratio still holds Illumina short reads ~10x cheaper than PacBio.
RNA-seq (i.e. analysis of a gene expression) is not possible with PacBio due to preferential sequencing of smaller fragments; shorter genes would always be shown to be more expressed. To be clear, it's possible to sequence RNA with PacBio (the keyword is iso-seq), but the analysis of gene expression is problematic.
It's way easier to extract fragmented DNA (concerns small non-model organisms; although recently a single mosquito was sequenced, so we can expect a further improvement)
other sequencing techniques as RAD-seq that allow genotyping with very little effort and cost, I have never seen anybody even considering using long reads for such genotyping
Genome profiling (assembly-less genome evaluation) based on kmer spectra analysis is not possible with PacBio data due to higher error rates. Conflict of interest: I am a developer of one of the tools for genome profiling (smudgeplot)
DNA in Formalin-Fixed Paraffin-Embedded samples is fragmented anyway (~300bp), therefore, there is no point in sequencing them with more expensive long-read technology (contributed by @mRoten)
I bet there will be a plenty other applications I am not aware of


Get residue position from Pymol Selection I can make a selection in Pymol using the resi parameter. Is there a way to get this resi data back out of a selection object? The only thing I found so far is how to get the position, as in the x/y/z coordinates from a selection.

For context, I have a table with a set of data points for each residue in the amino acid sequence of the displayed structure and I want to write a function where I can print the data based on the current selection. So my thought would be something like the following, where the content of the 'get_residue_positions_for_selection' function would be what I'm asking about.

def show_selection_position_data():
    residue_data = {
        "header": ['variation', 'rmsf']
        , "1": ['high', 0.1]
        , ...
    }
    selections = cmd.get_names('selection')

    if sele not in selections:
        print("No selection")
        return

    residues = get_residue_positions_for_selection("sele")

    for residue in residue_selections:
        print(residue_data[residue])

    return

Ok, after some more searching I found the answer.

You have to use the cmd.iterate function to iterate over the default selection ("sele", first parameter). With the second parameter you can push the variables you need into an array for later access. In this case I wanted the residue number so I used the variable 'resv'. A full list of the variables exposed by the iterate function can be found on the pymolwiki

Note that you can use a function call as the second parameter and if you need a more complex selection than 'sele' you can use selector.process in the iterate function.

Here is my final working function. Note that I used a set for storing the residue ids, since the iterate function technically iterates over the atoms and so the same residue will likely appear multiple times in your selection.

def show_selection_position_data():
    residue_data = {
            "header": ['variation', 'rmsf']
            , 1: ['high', 0.1]
            , 2: ['high', 0.2]
            , 3: ['low', 0.1]
            , 4: ['high', 0.05]
            , 5: ['middle', 0.2]
            , 6: ['high', 0.1]
            , 7: ['high', 0.1]
        }
    selections = cmd.get_names('selections')

    if "sele" not in selections:
        print("No selection")
        return

    userSelection = "sele"

    stored.residues = set()
    cmd.iterate(selector.process(userSelection), 'stored.residues.add(resv)')

    for residue in stored.residues:
        if residue in residue_data:
            print(residue_data[residue])
        else:
            print("{}: no data".format(residue))

    return


The effects of incomplete bisulfite conversion upon mapping efficiency This question has also been posted on Biostars

I have sequenced numerous multiplexed pools of BS amplicon-seq libraries derived from human samples on a MiSeq over the past few weeks. I have been utilising trim-galore and Bismark for alignment and am finding the mapping efficiency to be really low for two pools (55% & 30% respectively) both of which had a cytosine per base sequence of around 10-20% throughout the entire read in their fastqc files.

The high C content visible in the fastqc files makes me think the poor mapping efficiency is due to poor bisulfite conversion, as this would be expected to be close to zero is bisulfite conversion had actually took place.

I am new to the field so any help would be greatly appreciated.

Bisulfite conversion efficiency has no effect on the mapping rate in bismark and similar tools. The reason is that the reads are fully bisulfite converted in silico before alignment to minimize mapping bias.

I would suggest that you play around with the settings handed to bowtie2, such as using local alignment and modifying the --score-min option to allow more mismatches. Alternatively, you might try a different aligner like bwameth, which will always do local alignment.


Feature extraction methods that can handle inconsistent numbers of atoms for molecular dynamics I want to compare the protein dynamics
) pH 7 versus pH 3, or
) wild type versus mutant

The protein will have slightly different number of atoms at each condition, due to protonation or mutation, which is bad for machine learning. The alternative is to only consider the alpha-carbon atoms, but this will lose significant amount of information (e.g. side chain, salt bridge).

Is there a way in molecular dynamics to extract features from ensembles of different number of atoms?

Note: I have used common methods to decipher the difference between two conditions, including RMSD, RMSF, Rg, secondary structure, salt bridge, SASA, etc. They can explain the difference to some extent. However, I feel the selection of properties is "subjective", and there must be some hidden features not known. So I am seeking a machine learning method, that can "automatically" capture the most essential properties based on the coordinates as inputs.

This answer is no longer on topic. The question is about what is the best machine learning algorithm to use to analyse atomic coordinates from an MD trajectory in order to infer a novel properties.

  Also, it should be preface for any confused reader that unfortunately "feature" is a technical term in both genetics and statistics. For the former a feature is a region of DNA, RNA or protein that has a given property, such as open reading frame or PFAM domain. While in the latter a feature is a preferably independent property that forms part of a feature vector for different types of machine learning.

Heavy atoms are common

Assuming the pH difference only results in protonations, the protons differ between the two conditions, but not the heavy atoms. These will retain all the side chain information.

Merit of backbone only

(Section no longer relevant due to OP edit)

The common metric used to compare two structures is root mean square deviation (RMSD) of the alpha carbons to assess how much does the backbone differ and with analysis of MD simulations it is used to assess whether one is more dynamic than the other (e.g. this Bioinformatics SE question). The reason why the RMSD is done with the alpha carbons and not the minimum common subset of atoms is that backbone changes are large scale, while the side chain changes only affect locally —so their inclusion would needlessly dilute the score. The C prime and backbone nitrogens behave like the C alpha so there is no need for these.

Extraction of coordinates

Assuming by features you mean atomic coordinates, there are loads of ways. Using Pymol module for Python, one could do something like:

import pymol2
with pymol2.PyMOL() as pymol: #can be parallelised
    pymol.cmd.load('trajectory.gro', 'mytraj')
    pymol.cmd.remove('h.') #remove hydrogens
    pymol.cmd.intra_fit("name CA", 1) #align all steps in the trajectory
    for i in range(1, pymol.cmd.count_states('mytraj')):
        for atom in pymol.cmd.get_model(f'state {i}').atom:
            atom.coord #List[int int int]

One could change pymol.cmd.get_model(f'state {i}') to select state {i} and name CA for Cαs only.

Edit: Mutated residues

In the case of mutated residues or residue differences due to the two proteins being close homologues, "mutating" the differing residues to alanine, i.e. by stripping all heavy atom beyond the Cβ (keeping C′, Cα, H, N, O & Cβ) in the MD trajectory —or glycine if one is glycine. If they are homologues with different lengths, this approach would mean skipping residues that match up with gaps in the MSA. Even though not uncommon for statistically analyses, it does pose some issues statistically. It may unbury hidden residues, but length variant generally happens with loops.

In the snippet above, the pymol mutagenesis wizard could be used, but that is overkill, does not work in parallel and the approach does not port to other PDB reading modules in whatever language one wants. So say you have a list to_be_trimmed of integers of PDB residue indices (not pose residue indices):

if atom.resi not in to_be_trimmed or atom.name.strip() in ('C', 'CA', 'H', 'N', 'O', 'CB'):
      do whatever with atom.coord

In acidic environments, you get some funky things happening like exposed aspartate-glycine residue pairs forming isoaspartate or exposed aspartate-proline resulting in a cleavage, etc. For these cases, which are nearly never modelled, ignoring these values is the best option.

Edit: small comment on SASA

The investigated property list is very extensive and some reader may adopt them. So I thought I should mention that per residue SASA is generally converted to relative surface area, which is residue independent, by dividing it by the maximum SASA for that residue (values from Tien et al. 2013: 'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0, 'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0, 'W': 264.0, 'Y': 255.0, 'V': 165.0}), where RSA > 0.2 is an exposed residue.


Annotating conserved protein domains in bacterial assemblies I understand the basic principles under which RPS-BLAST operates: search a protein (rpsblast) or nucleotide (rpstblastn) query against a database of protein profiles, a "reverse PSI-BLAST" search. My question has to do with the practical application of this tool to annotate conserved protein domains in bacterial assemblies.

Let's assume I have downloaded NCBI's Conserved Domain Database (CDD) and would like to annotate conserved protein domains in a bacterial genome I have just assembled from Illumina reads. I see two potential approaches.

Run a gene prediction tool on the assembly, translate the predicted gene models, use rpsblast to search the translated sequences against the CDD profiles, and then transform the matches from local protein coordinates to global genome coordinates.
Use rpstblastn to search the genome sequences against the profiles directly.

The latter option would be more compute intensive, since it requires computing a 6-frame translation of each genome sequence. But perhaps it is more sensitive in the presence of small misassemblies?

In practice, is either of these approaches preferred over the other? Are there any considerations I'm missing?

The first option. Using a pre-existing tool, such as PROKKA etc.

Diversity. Bacterial diversity is staggering and species such as E. coli and Salmonella, which are closely related diverged when mammals and dinosaurs/birds diverged, so searching nucleotide sequences is a bad idea as you will get nothing. Especially as you are hoping to find distant sequences that may be missed by the commonly used annotation tools.

Operons. Bacteria, especially thermophiles, often group their protein functionally in operons. Therefore, you can use the principle of guilt by association to guess what a mysterious protein does. It is extremely powerful. Say you have Voet & Voet handy, a pathway with a hole and a mysterious protein you can guess easily what it does.

Note, if you have a protein that you are particularly interested in, I would suggest making your own PSI PSSM file. The automated system cannot read the literature or do some quick checks to see when to stop, so may be overly conservative with the clusters.


What tools can I use for a bacterial core/pan genome pipeline? I want to perform a genome comparison on a group of isolates. I want to look into two broad groups of taxa and compare the accessory genome in each group. I have been using prokka (v1.12) and roary (v3.8.2) to do this but it appears the accessory_binary_genes.fa file is actually an untrue representation.

Note: gene_presence_absence.Rtab does contain all the full presence/ absence for accessory gene sets. Despite this Im still unhappy with the nomenclature of the gene groups [issue for another day]

  [github issue 335] Your best to ignore the accessory_binary_genes.fa file. It is just for creating a quick and dirty tree with FastTree. The file itself is filtered to remove very common and not common variation to speedup the tree generation, hence the difference in numbers.

  The top 5% and bottom 5% are excluded. It is truncated at 4000 genes.

I've been looking into alternative pipelines and a new software BPGA looks promising. Does anyone have experience with this?

I essentially want a tool which will give me the core and accessory gene sets, without the noise from partial gene hits.

LS-BSR should be able to give you what you're looking for. See the article.

After you run the primary analysis, there is a simple documented workflow for splitting the pangenome into core and accessory, based on a user defined threshold. I'm the developer, so can help if you run into any problems.


Filter Trinity transcriptome based on RNASeq reads I have recently generated a genome-guided transcriptome with Trinity, and would like to apply an additional filter to exclude transcripts that don't have good support from the RNASeq reads. This is with the goal of trying to reduce the initial dataset down to something a bit more manageable (I have about 300k transcripts covering 250Mb in total, but would prefer about 1/10 of that number).

I see that Trinity has a workflow using Bowtie2 for evaluating the quality of assembled transcripts, but Bowtie2 is getting a bit old in the mapping world now. Are any of the super-fast transcript mappers (e.g. Salmon or Kallisto) appropriate to use for transcript filtering?

I used a combination of BUSCO and Salmon to filter transcripts based on their abundance in the RNASeq read dataset. The approach was roughly as follows:

Run BUSCO in short/transcript mode on the Trinity-generated sequences
Run Salmon on the Trinity-generated sequences, using the RNASeq reads that were used to generate the transcriptome
Merge the BUSCO full results table with the Salmon results table to determine distribution of the number of reads mapped for BUSCO genes
Choose an appropriate "real" read count threshold based on the read count distribution for BUSCO genes
Sub-select the transcripts based on this threshold, using the Salmon results table

Step 1 -- BUSCO

python ~/install/busco/BUSCO.py -i ../Trinity-GG.fasta -o BUSCO_Trinity-GG_nematodes -l ~/install/busco/nematoda_odb9 -m tran -c 10

Step 2 -- Salmon

~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon index -t Trinity-GG.fasta -i Trinity-GG.fasta.sai
~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon quant -i Trinity-GG.fasta.sai -1 reads-all_R1_trimmed.fastq.gz -2 reads-all_R2_trimmed.fastq.gz -p 10 -o quant/reads-all_quant -l A

Step 3 -- Results merge (using R)

#!/usr/bin/Rscript
busco.table <-
    read.delim("BUSCO/run_BUSCO_Trinity-GG_nematodes/full_table_BUSCO_T-BNOCFED_nematodes.tsv", comment.char="#",
               header=FALSE, stringsAsFactors = FALSE,
               col.names=c("Busco id","Status","Sequence","Score","Length"));
busco.table <- subset(busco.table, Status != "Missing");

count.table <- read.delim("quant/reads-all_quant_ISR/quant.sf",
                          stringsAsFactors=FALSE);

busco.expr.df <- merge(busco.table, count.table,
                       by.x="Sequence", by.y="Name", all.x=TRUE);

Step 4 -- Count thresholding (visual exploration)

pdf("busco_NumReads.pdf", paper="a4r", width=11, height=8);
options(scipen=10);
busco.numReads <- sort(tapply(busco.expr.df$NumReads, busco.expr.df$Busco.id,
                              function(x){
                                  exp(mean(log(x[x>0])))}), decreasing=TRUE);
par(mfrow=c(1,2));
plot(x=seq(0,1,length.out=length(busco.numReads)),
     busco.numReads, log="y", ylab="Number of Mapped Reads",
     yaxt="n",
     xlab="Proportion of Complete BUSCO sequences");
axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
rect(ybottom=10^(par("usr")[3]), ytop=10^(par("usr")[4]),
     xleft=0.9, xright=par("usr")[2],
     col="#80808020", border=NA);
plot(x=seq(0,1,length.out=length(busco.numReads)), xlim=c(0.9,1),
     busco.numReads, log="y", ylab="Number of Mapped Reads", yaxt="n",
     xlab="Proportion of Complete BUSCO sequences");
axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
invisible(dev.off());

Looks like a value of 50 will be suitable.

Step 5 -- Sub-select the transcripts

cat(file="HighCount_Transcripts_Num50.txt",
    subset(count.table, NumReads >= 50)$Name, sep="\n");

Then filter using a suitable FASTA selection utility:

pv Trinity-GG.fasta | ~/scripts/fastx-fetch.pl -i HighCount_Transcripts_Num50.txt > HighCount50_Trinity-GG.fasta


How to use SBOL (Synthetic Biology Open Language)? One of my students is working on the development of a SynBio tool to design and represent genetic circuits in a standardized format. I know there are some libraries out there.
What I am interested in finding out is to have a hands-on tutorial how to use an SBOL library (say, pySBOL) to represent genetic circuits in SBOL data and visual representations.

For pySBOL there is a tutorial you can follow in the documentation. Additionally, there is a youtube video walking through a similar tutorial, note the tutorial starts at about 1:09:00.
There are lots of additional tools and libraries out there such as SBOL Canvas (for visualisation), libsbolj (a java library for SBOL), and other libraries and tools listed on the SBOL website.


Linking of MeSH terms and genes Is there an alternative to mesh2gene for linking MeSH terms to genes?

The http://www.ncibi.org/gene2mesh.html site appears to be down for months now (click on "Launch Gene2MeSH via gene2mesh.ncibi.org") and I've been unable to find an alternative so far...

I think you could use the MESH.db package in Bioconductor to annotate your genes. I have never used, but you should be able to do something like:

library("MeSH.db")
select(MeSH.db, keys = myGenes, keytype = "ENTREZ", columns = "MESH")

It depends on the MeSHDbi package from which you could retrieve other data form that database if needed


Details of DESeq2 modeling a batch effect When correcting my data for a batch effect using removeBatchEffect, some of the gene expression values become negative.
When searching for differentially expressed genes, I do not use the data above, but rather model the batch using deseq2 (design=~Batch + Condition).
However, I started worrying. When DESeq2 introduces a batch into the model, does it allow for negative values "behind the scenes"?
If it does, I do not understand how that can make sense, in the context of RNA-Seq data.

DESeq2 uses the batch information (and everything else in the design) to produce offsets for its GLM. For a background on that please check how linear models work, e.g. using the StatQuest series of statistics videos over at YouTube.
It still operates on the raw counts. The same goes for the normalization factors.
removeBatchEffect fits a linear model to the data including the batch information and then subtracts the batch component from the counts (that is basically the baseline difference).
If you are interested in preserving the integer nature of counts and preserving zeros as smallest values after explicit batch correction you may want to check ComBat-seq() from the sva package. It operates on raw counts and returns batch-corrected raw counts which you could then normalize and calculate CPMs from (if you need batch corrected CPMs with no negative values). I find it useful and prefer it over removeBatchEffect as it avoids the unfortunate negative counts which sometimes messes up plotting scripts that expect zero as the smallest possible value.


DESeq2 - Error in lfcShrink: 'coef' should specify same coefficient as in results 'res' I have created a DESeq object (dds) and now I want to access the results and apply the lfc shrinkage on them.
res <- results(dds, alpha=0.05)
resultsNames(dds)
[1] "Intercept"                                                  "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages"
[3] "Cell_type_T_cells_vs_Alveolar_macrophages"

Here we see all the coefficients I got from the results function.
Now I try to use the lfcShrink() function on the second coefficient:
res_im_vs_am <- lfcShrink(dds, coef=2, res=res, type="apeglm")
Error in lfcShrink(dds, coef = 2, res = res, type = "apeglm") :
  'coef' should specify same coefficient as in results 'res'

it doesn't work, so I try with the name directly copy-pasted from the above code:
res_im_vs_am <- lfcShrink(dds, coef = "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages", res=res, type="apeglm")
Error in lfcShrink(dds, coef = "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages",  :
  'coef' should specify same coefficient as in results 'res'

Same error.
According to the ?lfcShrink for coef: the name or number of the coefficient (LFC) to shrink, consult resultsNames(dds) after running DESeq(dds).
So I do not understand why it is not working. Anyone can help?


Fastq-dump script download X spots or all Im trying to write a script where an optional input of -X flag can be used, or if that info is not available download all reads.

my script as follow:

# usage : bash dump_qc.sh <run id> <number of reads>

#run id variable.
SRR=$1

#maximum reads variable
X=$2

fastq-dump --split-3 -X ${X} ${SRR}

So, if I run the script

bash dump_qc.sh SRR1234 100

Will download 100 reads from the specified SRR file.

What I would like to do is to make the second argument optional, if not specified download the whole data set!

You don't need to copy your arguments into variables, you can just use $1 directly. Anyway, this should be enough:

#!/bin/bash
# usage : bash dump_qc.sh <run id> <number of reads>

if [ -z $1 ]; then
  echo "ERROR: at least one argument need to be given!" >&2
  exit 1
fi

if [ ! -z $2 ]; then
    x="-X '$2'"
fi

fastq-dump --split-3 "$x" "$1"


Is there a software or web server to calculate thermodynamic parameter enthalpy change of cellulase? I want to calculate thermodynamic paremeter enthalpy change of cellulase, but I don't know what software can do this work. Please tell me, thank you very much !

You could use PISA assuming you have a PDB file of cellulase? Sorry it's a short answer and should be a comment, I don't have the reputation, however.


data visualization RNAseq : scaling data for PCA and cluster dendogram I have count data from a RNAseq experiment (2 samples are from normal cells and 3 samples are cells with a disease), and the data is already standardized by trimmed mean of M values (TMM).
I want to do some plots: biplot of Principal Component Analysis (PCA) and a cluster dendrogram to see if the samples normal/disease are well separated (there is a clear difference between them).
Since the data is already standardized (TMM) should I scale and center the data prior to perform PCA and cluster dendrogram??
thank you!

Scaling (or centering) makes the genes comparable: Putting the expression levels of genes in the same scale (i.e. between 0 and 1) sustains that all of your genes contribute equally to the PCA or distance calculations. On the other hand, without this step, such calculations would be dominated by the highly/lowly expressed genes. I get "better" results with scaling.


RcppArmadillo installation error When I am trying to install the RcppArmadillo Bioconductor library in R, it is giving me the following error:

/usr/bin/ld: cannot find -lgfortran
collect2: error: ld returned 1 exit status
make: *** [RcppArmadillo.so] Error 1
ERROR: compilation failed for package ‘RcppArmadillo’
* removing ‘/home/nikita/R/x86_64-pc-linux-gnu-library/3.2/RcppArmadillo’

I am using Ubuntu 12.04, R 3.2.5 and g++-5.4. Any help would be greatly appreciated.

On Ubuntu 14.04 there is a direct route to installing RcppArmadillo. You could perhaps try the same syntax.

sudo apt-get update
sudo apt-get install r-cran-rcpparmadillo


why in 3' quantseq should we use reference genome and not transcriptome? I am studying about RNA seq. I have found that in 3'prime rna seq (quantseq, tagseq) we should use a reference  genome  with    good    annotation  (plus known UTRs). My question is why?
Why could not we use reference transcriptome?
Thank you!

It would be great if you can elaborate on where you "found" this information.
In the transcriptome annotation you have a defined number of transcript per genes, but there might be others that are present in your sample and not in the annotation. For example you have treated the cells with some extreme condition that changes the 3' processing machinery.
With 3' quantseq, you are left with the read fragment that is at the end of your transcript. If the ends of these transcripts are not part of the sequences in the annotated transcriptome, these will be unaligned.
These instances might be rare if your transcriptome annotation is good, but let's say if you find that a gene of interest which should be expressed is not, and you suspect it might be due to a different isoform used, you can always check it with the genome alignment (see where reads might be falling in the gene). With a transcriptome alignment, this will not be possible.


cDNA and alignment mapping I am confused with RNA seq alignment. My understanding is that after Mature mRNA is isolated from the cell, it is then fragmented and using reverse transcriptase enzyme a cDNA copy is created which is further amplified. This cDNA is then mapped with the reference genome. If the cDNA does not have introns how does it map with the reference genome.

Can someone please help

tl;dr

Spliced aligners just make many small matches along the read pair or single read to determine splice junctions and output an alignment for the whole read, as long as it matches somewhere in the genome.

mRNA and translation

In the cell, after transcription, introns are spliced out of RNAs to become mature RNA. The protein-coding sequence does have sequence that is chimeric from cis-located portions of the same chromosome that do not actually exist.

              |<<exon 1>>|  intron    |<<exon 2>|
Genome: ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

mRNA  : ----ABCDEFGHIJKLYZ123456789aaaaaaaaaaa

As you said, if we only even consider small stretches of DNA, the sequence JKLYZ does not exist in the genome, only in the mRNA. So, we need a special way to map mRNA-seq reads to the genome that correctly splits at splice junctions.

How do we map mRNA-seq reads to a genome?

As a result, parts of the mRNA-seq reads will align well to one part of the genome, but only for part of the read. This requires a different alignment algorithm that considers multiple local high-score alignments for different parts of the read.

For example, see how a genomic aligner versus a splice aligner might treat the same mRNA-seq read on our example genome.

                         |<<exon 1>>|  intron    |<<exon 2>|
Genome:            ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

spliced aligner  : ------ABCDEFGHIJKL------------YZ123456789
genomic aligner  : ------ABCDEFGHIJKL

There is a possibility that the genomic aligner will not recognize that the second half of the read is an exon due to the large intron gap, and it may clip the alignment. Spliced aligners just take this into account and output several alignments for the same read regardless of gap size. The STAR algorithm, for example, uses many small matching pieces of a single read to determine where the splice junctions are. HiSat2 uses small approximate matches using a FM-index to find seeds for spliced alignment.


How can I compare structures of proteins responsible for dehydration tolerance in plants? For my undergraduate research project, I'm trying to compare the structures of some proteins, specifically transcription factors crucial for dehydration tolerance in plants, in the hope of finding differences that make each protein more or less dehydration tolerant than the others. i want to compare the selected transcription factors at their structural level. structure of one transcription factor from different plant species such as rice,sorghum,maize etc. are considered. ultimately i'm hoping to find a relationship between these structures of proteins and their respective level of dehydration tolerance in each plant.

Is there any database or a site where I can find information regarding this?

Have you explored in PDB?. If you already have a starting sequence, you can use it as input to search with phmmer using PDB as database. Once you have a curated set to analyse, you can visualize and align the structures with VMD, for example.


Merging sequencing data for ChIP-seq experiments I need to merge sequencing data from different sequencing runs but for the same ChiP-seq library (HiSeq 2000).

Are there any potential advantages or disadvantages when merging files at .fastq or .BAM stage (alignment with Bowite/1.1.2)?

I don’t think it matters. Both are easy to merge (BAM via samtools merge, and (gzipped) FASTQ via cat), and neither method has specific disadvantages, unless your FASTQ files are sorted for some reason (but they generally shouldn’t be).

One advantage of keeping the FASTQ files separate is that it makes it slightly easier to parallelise the mapping step: just run the mapper in parallel on the separate FASTQ files. Although bowtie has an option (-p) for this, throughput from that is slightly worse than running the mapping on split files.


Aligned base strand in pysam `pileupcolumn` Is it possible that on the same pileupcolumn object (i.e a specific aligned base) I find the same base but aligned to different strands? For example, can I find a coverage of 50 'A's that align to strand '+', and on the same time 50 'A's that align to strand '-' within the same index? I expect that if 'A' is aligned to strand '+', I will find 'T' aligned to strand '-'.

As @terdon mentions, it's certainly possible to have an A aligned on both strands. If a read has a mismatch at the location of interest and that mismatch is the complement of the reference base, then you are likely to see A in both strands. However, it is very unlikely you will see the same base at high coverage in both strands. For example, if your total coverage is 55, and you see 50 As on the + strand and 50 As on the - strand, it's likely that the pileup object is reporting the reverse complement sequence of the reads that align to the - strand for your convenience.


Differentially expressed genes analysis in Seurat For the differentially expressed genes analysis, is it possible to check for DEGs based on the levels already identified in the object? For example, my dataset contains cells from 11 subjects, which are marked as so (orig.ident). I want to analyze the differentially expressed markers between them. Will it be accurate to do so? Or should I follow the clustering analysis in Seurat? I want to see the difference in genes between the cells according to the subjects but I'm afraid that may not be the right way to do so.
  head(BC@meta.data)
         orig.ident nCount_RNA nFeature_RNA percent.mito
 BC01_02       BC01   998944.0         9189 0.2387108038
 BC01_03       BC01   999696.4         9548 0.2925345423
 BC01_04       BC01   999057.5         7440 0.0009893724

Is clustering analysis absolutely necessary to see the DEGs between them? or can I directly move to FindMarkers based on the levels (orig.ident)?

I'm not quite sure what you are asking, but you can use anything in the metadata as a grouping with FindMarkers.
DefaultAssay(combined.all) <- "RNA"
Idents(combined.all) <- "orig.ident"

for (sample in unique(Idents(combined.all))) {
  clustermarkers <- FindMarkers(combined.all, ident.1 = sample)
  print(paste("Cluster", sample))
  print(head(clustermarkers, n = 10))
}