Replicate papers (beginner) I have a strong mathematical and computer background (also with machine learning) and would like to switch to bioinformatics. I know a fair amount of biology. But I find it hard to find real world practical training in this field. So I would like to try to replicate some bioinformatics papers (with their provided data). Good beginner papers for me would be shorter ones where the biological problem is well stated and the problem is also fitted in to the bigger scope (or is just easily understandable). Do you have any papers to recommend?

Bioinformatics is a big and very heterogeneous field with a lot of variation, so it's hard to recommend something that covers everything. That said, I think the best you'll do is to build familiarity with some of the common tools with some canned analyses.

One possible option would be some of the scRNA-seq workflows from the Pachter group, they seem to be fairly conscientious about documentation. Here is one such canned analysis, here is the paper, here is a blog post about it that suggests next steps and discusses at a high level.

scRNA-seq is just what's been trendy for a few years, next year it will be something different. But it at least gets you familiar with some of the file formats and concepts such as alignment, etc.

Hope that helps.


Use of Electronic Phenotype in EHR May I know what's the use of Electronic Phenotyping using EHR data?

I did refer this link but have few questions

I understand that Phenotype is a set of criteria that you apply on your EHR data to select patients of interest from EHR database.

However what I am trying to understand is the difference between ordinary SQL query and electronic Phenotype

If both are same,

1) what I would like to know is the use of electronic phenotyping? How it is used?

2) Why do we have pheKB repository and other repositories? What's the use in validating them in multiple sites?

3) Is it like identifying patients real time? I am not able to get the usefulness of this. I understand that we select patients of our interest (cohort), but i can just write an ordinary SQL and send it to my EHR database to identify patients. but what's special about Electronic Phenotype?

4) I also see that they use multiple conditions (Lab tests, Medications, Diagnosis etc) to identify patients. What am I missing here?

Can someone help me by giving an simple example? I am new to healthcare and trying to learn. Would really be helpful

I understand that Phenotype is a set of criteria that you apply on your EHR data to select patients of interest from EHR database.

No, a phenotype is a way of behavior or other observed characteristic of a person, resulting from a combination on its underlying biology and the environment it's in. Think things like, height, running speed, strength, etc.

In the context in which you're interested, this would be more properly, "the clinical presentation of a disease or disorder".

Step one of diagnosing a patient is defining exactly how they present (i.e., exactly what their maladies, aberrant electrolyte levels, genetic mutations, etc. are).
An algorithm that only works on the test population is useless. Google "model overfitting".
The point of an algorithm is to take the EHR data and progress to a diagnosis. This isn't real time, you'd just run a batch of algorithms as test results come in.
You need patients in your test/validation sets to have the required inputs for your algorithm.

FYI, I don't know of any regulars on this site that deal with EHR data. If there's a community of such people on the web they must be elsewhere.


Download sequences of isoforms I want to collect all isoforms of all genes (Fasta/Fastq file of nucleotide and protein sequences) Arabidopsis Col-0. I am wondering if there is a straightforward way to download the file from any database.

I am new to sequence analysis, so I apologize if asking a very basic question.

Col-0 is, to my knowledge, the standard reference Arabidopsis thaliana ecotype, so the easiest way to get the sequences for all genes is to just download them from Ensembl. The nucleotide sequences are available here and the protein sequences here.


How to sort multiple FASTA files based on their content? I have around 10,000 FASTA files of Influenza A virus.
These files contains sequences of each of the 8 segments of the viral genome and I want to separate these files into different locations based on the content of these FASTA files.
In each FASTA file for each segment,the first line has the segment number, for example.

KM368312.1 Influenza A virus (A/swine/Shandong/01/2009(H1N1)) segment 3 polymerase PA (PA) and PA-X protein (PA-X) genes, complete cds

To be clear I want

file1.txt has segment 1
file2.txt has segment 2
etc ...

I want to ensure all segment 1 sequences are grouped into one folder, and each file is grouped according to its geographic origin. The geographic groupings are mirrored for all 8 segements and each placed into their own directory.

The general idea is:

iterate over all files

read the first line of each file, extract the segment bit
move the file based on that

In Python, this can be done in many ways. Here’s one way:
import os
import re

def find_segment(id):
    return re.search(r'segment (\d)', id).group(1)

for segment in range(8):
    os.mkdir(f'segment-{segment + 1}')

for entry in os.scandir():
    if not entry.is_file():
        continue

    with open(entry.path, 'r') as file:
        segment = find_segment(next(file))
    os.rename(entry.path, f'segment-{segment}/{entry.path}')

Note the absence of error handling — this code assumes that every operation succeeds; empty files, or files without “segment X” in their FASTA ID, will trigger an error, as will the usual IO failures (missing permissions, etc.).

But most people would probably use a shell scripting language instead of Python to solve this. For example, here’s something equivalent in Bash (untested):
mkdir segment-{1..8}

for file in *.fa *.fasta; do
    segment=$(sed 's/.*segment \([[:digit:]]\).*/\1/; q' "$file")
    if [[ -n $segment ]]; then
    mv "$file" "segment-$segment"
    done
done


Low percentage of reads with consistent barcodes in Split-Seq run I am following this answer to select reads with allowed barcodes from a Split-Seq run. In particular, I have a whitelist for each round of barcodes (3 rounds) and I filter out reads which have more than 1 Hamming distance from allowed barcodes.

However, after filtering, the number of reads remaining are only the 0.7% of the original number of reads before filtering (169k from 22M). Is this normal, or what may be the underlying issue here?

In my case the problem was due a wrong set of allowed round1 barcodes, therefore the filtering pipeline was removing most part of the reads. In particular, with the right barcode1 whitelist the total number of reads after filtering is 17M (about 100 times more).


FASTA and PDB: How to specify chain? For proteins that have multiple chains (e.g., 1EMS), is there an easy way to specify which chain I want to use for blastp?
I cannot imagine that I am the first person to have this problem, but so far my only solution has been to write my own script that gets the chain I want from the .pdb file and convert that to a .fasta file.

The NCBI BLAST can be used with the PDB DB (which NCBI has). The PDB codes are stored as 4 letter codes underscore chain, e.g. 1GFL_B. The catch is segment identifiers —but generally they are the same peptide so shouldn't be an issue.
You can search specifically the PDB DB in NCBI (not the RCSB PDB) by setting the database to PDB.


No counts for added gene in cellranger (scRNA-seq) I have a set of scRNA-seq samples enriched with FACS for cells expressing a specific gene reporter (TdTomato). In particular the gene I want to report has positive counts in the resulting matrix for 97% of cells.

I followed CellRanger documentation to create a custom reference. I added TdTomato on an additional chromosome and I added an entry in the gtf:

TdTomato    .   exon    1   1431    .   -   .   gene_id "TdTomato"; transcript_id "TdTomato"; gene_name "TdTomato"

The command I'm using is:

cellranger-2.1.1/cellranger mkref --genome=refdata-cellranger-mm10-2.1.0.tdtomato --fasta=genome.tdtomato.fa --genes=genes.tdtomato.gtf

where genome.tdtomato.fa is the original reference genome (mm10, as provided by cellranger) with an additional chromosome at the end with header:

>TdTomato

and genes.tdtomato.gtf the original reference genome gtf with the additional line indicated previously at the end of the file.

I actually tried both on the '+' and on the '-' strand, however I guess this should not affect the read mapping, since TdTomato is on an independent chromosome.

I ran cellranger count using the custom reference:

s="SAMPLE"
cellranger-2.1.1/cellranger count --id="$s"_tt --transcriptome=refdata-cellranger-mm10-2.1.0.tdtomato --fastqs=original_fastqs_directory/"$s"/fastqs --sample="$s"

I am using the default parameters of CellRanger.

However, I found that less than 1% of cells seems to express TdTomato. If I use the '+' strand I get slightly more counts (2.4% of cells have positive TdTomato).

What I checked:

With IGV I can see that there are reads mapping to TdTomato
Using samtools I am able to see the same number of reads mapped on TdTom either with the + or - strand
However, following this documentation I noticed that by executing the CellRanger filtering steps, more reads are filtered in the case of the - strand

This reduces the number of reads used for counts to 44 (2.4% of the original 1807 reads, - strand), 240 (13.3%, + strand) for TdTomato, 71620 (21.2%) for the gene reported by TdTomato. Interestingly, there are fewer reads for TdTom than counts, so there may be some kind of normalization afterwards?

The gene reported by TdTomato is longer, 11394 bp than TdTom, 1431 bp. The total length of the exons for the two genes is 5349 vs TdTom: 1431. Does this influence exponentially the amount of reads mapped to the two genes?

The first column of your GTF must correspond exactly to the name of the sequence in your FASTA file. Typically this is the chromosome so a GTF line:

chr1    hg38_knownCanonical exon    11106535    11262507    0.000000    .   .   gene_id "gene1"; transcript_id "tx1";

This must have a corresponding sequence in the FASTA file:

>chr1
ATCG ATCG ATCG ATCG ATCG ATCG...

Therefore you need to add an additional sequence to your FASTA file with >TdTomato. This must have exactly the same sequence as the construct that you are expecting reads for. The run cellranger mkref and cellranger count on these FASTA and GTF files. Cellranger also can take a GTF with “CDS” as a feature in addition to “exons”. The STAR alignment algorithm used by cellranger is splice aware so it can handle split reads but you will need to ensure the GTF corresponds to exons of the mRNA/cDNA rather than the full sequence.

10x Genomics is a platform with a high dropout rate, especially for lowly expressed genes. You should bear this in mind for what you expect to observe. It may be worth comparing the mapping rate of your construct to housekeeper genes to see if this is the reason for the low counts. Since whole mRNA molecules are captured on the beads before preparing the library, there should not be a difference in UMI counts between long and short genes.


HTSeq-count is returning 0 for every gene, instead of expression value (Ported from stack overflow)

I'm trying to summarize gene count using htseq-count; and it's returning 0 counts at every gene. I'm not sure what I'm doing wrong; I think it has to do with the gene flag I'm using.
I've tried using the GTF for Arabidopsis TAIR 10: which I got from: ftp://ftp.ensemblgenomes.org/pub/release-46/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gtf.gz as well as the GFF file available on Araport.

I'm using various BAM files I know to have gene expression values.

I've tried both locally and on galaxy
The local commands I've tried so far are:

python3 -m HTSeq.scripts.count --idattr=gene_id #bam #gtf > summary.txt
python3 -m HTSeq.scripts.count -i gene_id #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i locus_tag #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i ID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i SeqID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i Parent #bam #gff > summary.txt

I've done the same on galaxy; inputing the -i value in.

Any help as to where I'm going wrong would be great.

Thanks!

Maybe the contig names in your BAM files don't match those in the GFF file? Ensembl uses contig names like 1, some other reference sources use chr1.

Generally, for ensembl GFF files, HTSeq's default value for idattr should be fine though.


Bulk RNA-Seq Read Length Normalization across different samples I have 20 samples out of which 14 are 100 bp in length and 6 are 150 bp. Is there a way to normalize the read length for cross-sample differential expression comparison? What would be the best way to approach this?
Note: The samples are NOT from reference/model species.

The simplest thing to do is to trim the 150 bp fasts so that they are 100 long.  I don't think there is an easy way to correct for the fact that the 150 bp long reads will have a higher unambiguous rate of alignment and gene assignment than the 100 bp long reads.
If you have a mix of all the experimental conditions across all the lengths, you can include length as an element in the design, so that DESeq2 may be able to model the difference between the two batches.


Can I automate CLUSTALO and output alignment sequence identity? I've detected homology between targets of ligands in drugbank and proteins in the proteome of a pathogen. I've parsed the output very rudimentary and calculated my query coverage. This exists in an excel file (see attached)

I'd like to further filter these hits. For query coverages > 70% alone I have 1500+ hits. Some redundancy exists and some hits are for drugbank targets for which no data has been deposited in the PDB.

So there are ways of further filtering these hits, but I'd like to do pairwise alignments first before I do anything else, as I think this is the most rational filter I can think of.

I do not want to filter based on lack of PDB data necessarily, as when I check my high query coverage hits, I'm getting predictions that are actually backed up by in vitro data. This includes high coverage hits lacking x-ray data

My question:

I have CLUSTALO installed on my desktop: How do programmatically send it jobs if I have the Drugbank sequences and the pathogen sequences in the same file? I would also like to parse out the pathogen IDs from the excel file. There’s 1500 and they are a bit scattered throughout the columns.

0.980634         -                    572         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-263             874.5          0.0            1        1      6.6e-267  1.4e-263   874.2      0.0       5         562    15     570    12     572    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005678  PE=3             SV=1

0.950704         -                    580         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-258             858.0          0.0            1        1      6.8e-262  1.5e-258   857.7      0.0       6         546    10     560    7      578    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005014  PE=3             SV=1

0.950704         -                    691         sp|P06672|PDC_ZYMMO                         -                     568                               8.2e-162             538.3          0.1            1        1      2.4e-164  5.1e-161   535.6      0.1       1         541    89     650    89     658    0.89   Acetolactate                            synthase                  OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_006830  PE=3             SV=1

0.492958         -                    585         sp|P06672|PDC_ZYMMO                         -                     568                               8e-87                290.7          0.0            1        2      3.3e-83   7e-80      267.8      0.0       3         283    2      287    1      292    0.94   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_001897  PE=3             SV=1

0.107394         -                    585         sp|P06672|PDC_ZYMMO                         -                     568                               8e-87                290.7          0.0            2        2      1.7e-08   3.6e-05    21.1       0.0       413       474    436    505    376    583    0.72   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_001897  PE=3             SV=1

So, as the first reply said, there's sort of multiple questions here.

For filtering things out from a table, I would absolutely use pandas. https://swcarpentry.github.io/python-novice-gapminder/08-data-frames/ Pandas is a way of managing dataframes in python. It looks like what you want from your table is to remove certain rows, if the content of those rows at a certain column matches certain IDs from a list (pathogen IDs). That's the sort of thing the pandas library is good at.

For managing workflows, I use either bash or snakemake. Snakemake is a python thing. It has some advantages, such as being able to start from an intermediary file (for example, if a workflow is "sequences -> multiple alignment -> phylogenetic tree", it can easily start from the multiple alignment). Also deals with branched workflows quite well. Can be easily installed with conda (anaconda is absolutely a lifesaver for bioinformatics work imo). I still use bash for simple workflows, because it's slightly simpler/faster.

For potential filtering - If you want to keep the top hits as they are, but filter the lower-scoring hits based on whether they have a PDB ID or not, I'd personally look into Entrez queries. They're a way of making your search more specific. Here they're shown for WebBlast, but they work quite well for command line blast too. https://www.youtube.com/watch?v=bxx5uaKjMa8


Reducing a string to a certain amount of the first characters in R So I have a list of several thousand RNA transcripts, and I am trying to find out how they are spliced. Right now I am running matchLRpatterns() from the Biostrings package with a max gap length of 0, after running a matchpattern function to categorize the transcripts by donor sites (where the first cut in an RNA transcript is made to cut out introns). But applying that to several thousand transcripts is quite time consuming, when you have 5 donor sequences to match to 9 acceptor sequences in 100 nucleotide long transcripts. Instead, I am thinking to sort by donor, remove all that comes before the splice junction by converting it to a string then using gsub, then using matchpattern again, but I only want to check the first 8 nucleotides after the splice junction for the consensus sequence that follows my list of acceptor sites. If I just run matchpattern without shortening the length of the transcript, I will just find acceptor sites in the RNA transcript, not the acceptor site at the splice junction.

I was thinking lapply the gsub function to the transcripts to remove up to the splice junction, then lapply the strsplit function to make a list of lists of groups of 1 character strings, lapply the head function to get the first 8 single character strings in each list, lapply the paste function to make them one string again, lapply the DNAString function to make the sequences readable by matchpattern, then lapply the matchpattern function to check the acceptor sequence at the splice junction. Is there a better way? This feels a bit clunky and roundabout way of doing it. I know gsub can remove characters after a certain character or sequence of characters, but I am not sure what the first sequence in each transcript is at the splice junction. Can I still use gsub, but tell it to remove characters after a certain number of characters? Or better yet, is there a way to get matchpattern to only check the first x nucleotides in a sequence?

Is there any particular reason why you want to do this via pattern matching in R? This sounds like something that SuperTranscripts and Lace might be useful for:

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1284-1

  Finally, Lace will annotate each superTranscript with blocks either by using the graph structure itself or alternatively using a script called Mobius (packaged with Lace) which infers the annotation from spliced junctions discovered when mapping reads to the superTranscript.


Can I get longer contigs by changing MEGAHIT settings? My current settings are megahit -r input.fq --num-cpu-threads 32 --min-contig-len 300 --presets meta-large -o output. I picked the 'presets meta-large' because I am dealing with a quite large and diverse metagenome out of which I only need some genes.

I've got a contig of interest ~19kb, with a depth of 16. Could I make this contig larger with more relaxed settings? Help would be appreciated!

You can use the option --k-step 10 and --kmin-1pass for low-coverage datasets.


Tumor sample and heterozygous SNVs: what is the genotype of the normal cells? In the following picture, you can see that in normal cells (yellow) we have two copies of each gene, while in tumor cells (violet), we just have one copy due to a clonal mono-allelic deletion. I do not really care about the genotype of the two genes, I just need to know that they are heterozygous.

On the other hand, suppose we are talking about somatic heterozygous SNVs and therefore we do not consider deletions anymore. In this case I would like to know the genotype (AA, AB or BB) of tumor and normal cells. Based on my understanding, I expect the AF of the SNVs in the tumor cells to be roughly 0.5 (since I ask for heterozygous SNVs). But what about normal cells? What is their genotype? If I ask for heterozygous SNVs, does that imply that normal cells are homozygous for that alleles?

Your initial information/image and actual question seem pretty unrelated to me, but I'll do what I can.

If you have a heterozygous SNV in a malignant population of cells, it really depends on when the mutation occurred for you to know the genotype of wildtype cells. It could have been there since birth (a germline mutation passed down from the parents) or have been acquired throughout the person's life and clonally expanded (somatic mutation).

Unfortunately, you have no way to differentiate between these two types of SNVs just by sequencing the tumor cells. If you want to know true somatic mutations that are found only in the tumor cells, you need to sequence a matched normal control from the same patient (preferably from the same cell type). This is why most large cancer sequencing studies include not only tumor samples, but normal cells from the same patient - see TCGA or ICGC for examples.


Building STAR Genome Index for nanopore RNA sequencing I am aligning a dataset of 1,000,000 reads oh human mRNA sequenced on Oxford Nanopore Technologies' MinION, and would like to use the STAR aligner, using the parameters recommended by Pacific Biosciences for long reads.

According to this Google Groups thread, in setting up the genome index for short reads, the parameter sjdbOverhang should be set to 1 less than the read length. Obviously, with long (mean 1.7Kb, max >50Kb) reads, this doesn't make sense.

Running STARlong --runMode genomeGenerate without setting sjdbOverhang sets the parameter to 0. Does anyone know what this means, how it might affect mapping, and what I should set it to for my long reads?

The parameter is used to determine how much sequence STAR indexes on each side of a splice junction to improve its alignment accuracy.  For very long reads, this may not be ideal. I am not sure if STAR is capable of including multiple splice junctions since a long read is more than likely to span more than one.

It may be worthwhile to consider aligning your reads directly to transcripts and taking the reads that align poorly and using them to look for novel transcripts.


Database of position weight matrices for protein motifs? I am trying to identify proteins that carry a consensus target sequence for a kinase. Usually if I am working with, for e.g. transcription factors binding to DNA I would use position weight matrices (PWMs) from one of a number of databases of PWMs to search the sequences with something like FIMO.

However, I don't know of any databases that contain protein motif PWMs. Specially I am looking for one for yeast Cdc28, but it seems like it would be a generally useful resource that I can't believe doesn't exist. For example, while the MEME suite was originally developed to identify protein motifs, and ouputs PWMs, all of its databases refer to DNA or RNA.

PFam has HMM models for protein domains, this is almost certainly what you're after. You can download all PFam-A models from ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz.

Then, you're looking for the "Pkinase", PF00069 family (I found this simply by searching the PFam website for "kinase"), so extract it from the file:

zcat Pfam-A.hmm.gz | awk -vRS='//' -vORS='//' '/PF00069/' > PF00069.pfam

Then, format it for HMMER format (you will need to install HMMER)

hmmpress PF00069.hmm

And now you can use it to search your sequences (fasta or flat file, probably others; see the documentation for details):

hmmscan PF00069.hmm prots.pep

I tried this on a fasta file with the sequences of CBK1_YEAST, CSK21_YEAST (kinases) and FLNA_HUMAN (non-kinase control) and got:

# hmmscan :: search sequence(s) against a profile database
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query sequence file:             prots.pep
# target HMM database:             PF00069.hmm
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Query:       sp|P53894|CBK1_YEAST  [L=756]
Description: Serine/threonine-protein kinase CBK1 OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CBK1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    2.7e-64  203.2   9.0    3.5e-46  143.8   0.2    3.1  3  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 ?   -2.9   3.2      0.19      0.19     198     226 ..     222     243 ..     191     277 .. 0.57
   2 !  143.8   0.2   3.5e-46   3.5e-46       1     153 [.     352     505 ..     352     522 .. 0.95
   3 !   63.6   0.0   9.8e-22   9.8e-22     154     264 .]     568     672 ..     551     672 .. 0.83

  Alignments for each domain:
  == domain 1  score: -2.9 bits;  conditional E-value: 0.19
                           TTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
               Pkinase 198 sgekgkekvekeldqlekilkilgetkek 226
                           + +++++++++      + ++i ++++++
  sp|P53894|CBK1_YEAST 222 Q-QQQQQQQQQ------QHMQIQQQQQQQ 243
                           1.111111111......122222221222 PP

  == domain 2  score: 143.8 bits;  conditional E-value: 3.5e-46
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHH.HHHHHHHHHHHHHSSTTB--EEEEEEESSEEEEEEE--TTEBHHHH CS
               Pkinase   1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkke.kkvlrEikilkklkhpnivrlyevfeekdelylvleyveggeladl 88
                           +++++ +G+G+fG+V  + +k+tgki+A+K++ k+++ kk++  +v +E  +l+ +++p +v ly++f++ ++lyl++e+++gg+l+++
  sp|P53894|CBK1_YEAST 352 FHTVKVIGKGAFGEVRLVQKKDTGKIYAMKTLLKSEMYKKDQlAHVKAERDVLAGSDSPWVVSLYYSFQDAQYLYLIMEFLPGGDLMTM 440
                           67899********************************999999********************************************** PP

                           HHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEETTTEEEE-SGTTSEECCSSCS CS
               Pkinase  89 lsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLidekgelKitDFGlakelesssk 153
                           l r + ++e+ +++++++ + ++e +H+ g+iHrD+Kp+NiLid +g++K++DFGl++ ++++++
  sp|P53894|CBK1_YEAST 441 LIRWQLFTEDVTRFYMAECILAIETIHKLGFIHRDIKPDNILIDIRGHIKLSDFGLSTGFHKTHD 505
                           **999*****************************************************9999876 PP

  == domain 3  score: 63.6 bits;  conditional E-value: 9.8e-22
                           BCTCCSCGGGS-HHHHTTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH......TTS-HHHHHHHH CS
               Pkinase 154 ltsfvgtreYlAPEvlkeneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkeklpeaselseeakdllk 242
                           + s vgt++Y+APE+  +++y++++D+WslG i+ye l g ppf +e+ +e+++k+++          e+++++p + ++s ea dl++
  sp|P53894|CBK1_YEAST 568 AYSTVGTPDYIAPEIFLYQGYGQECDWWSLGAIMYECLIGWPPFCSETPQETYRKIMNF---------EQTLQFPDDIHISYEAEDLIR 647
                           6789*****************************************66665533333333.........557799*************** PP

                           HHT-SSGCCSTT....HHHHHTSGGG CS
               Pkinase 243 kllkkdpkkRlt....aeellqhpyl 264
                           +ll+    +Rl+    a+e+++hp++
  sp|P53894|CBK1_YEAST 648 RLLTHA-DQRLGrhggADEIKSHPFF 672
                           ***875.678877779********98 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (756 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.00
# Mc/sec: inf
//
Query:       sp|P15790|CSK21_YEAST  [L=372]
Description: Casein kinase II subunit alpha OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CKA1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    7.5e-62  195.1   7.2    3.2e-55  173.4   0.7    2.1  2  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 !   25.5   1.1   4.3e-10   4.3e-10       1      60 [.      40      94 ..      40     115 .. 0.87
   2 !  173.4   0.7   3.2e-55   3.2e-55      56     264 .]     128     363 ..     101     363 .. 0.92

  Alignments for each domain:
  == domain 1  score: 25.5 bits;  conditional E-value: 4.3e-10
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHHHHHHHHHHHHHHHSSTTB CS
                Pkinase  1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkkekkvlrEikilkklkhpni 60
                           ye ++k+G+G++++V+++v+ +++ ++++K++k  k+kk k     rEikil  l+ +++
  sp|P15790|CSK21_YEAST 40 YEIENKVGRGKYSEVFQGVKLDSKVKIVIKMLKPVKKKKIK-----REIKILTDLSNEKV 94
                           78999****************************99888855.....7******9998877 PP

  == domain 2  score: 173.4 bits;  conditional E-value: 3.2e-55
                            S.STTB--EEEEEEE..SSEEEEEEE--TTEBHHHHHHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEE.TTTEEEE CS
                Pkinase  56 k.hpnivrlyevfee..kdelylvleyveggeladllsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLid.ekgelKi 139
                            + h ni++l++++++  +++  lv+eyv++ ++  l+    +l++ e+++++ ++l++l+y+Hs+gi+HrD+Kp+N++id ++++l++
  sp|P15790|CSK21_YEAST 128 NgHANIIHLFDIIKDpiSKTPALVFEYVDNVDFRILYP---KLTDLEIRFYMFELLKALDYCHSMGIMHRDVKPHNVMIDhKNKKLRL 212
                            35*************8878888*********9988887...5**************************************66677*** PP

                            -SGTTSEECCSSCSBCTCCSCGGGS-HHHH.TTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
                Pkinase 140 tDFGlakelessskltsfvgtreYlAPEvl.keneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkek 226
                            +D+Gla+ ++ + +++ +v++r +  PE+l ++++y+ + D+Ws G +l++++++++pf   +g++    ++dql+ki+k+lg+++++
  sp|P15790|CSK21_YEAST 213 IDWGLAEFYHVNMEYNVRVASRFFKGPELLvDYRMYDYSLDLWSFGTMLASMIFKREPFF--HGTS----NTDQLVKIVKVLGTSDFE 294
                            ******************************999***************************..3443....789**********77766 PP

                            ....TTS...............................-HHHHHHHHHHT-SSGCCSTTHHHHHTSGGG CS
                Pkinase 227 lpeasel...............................seeakdllkkllkkdpkkRltaeellqhpyl 264
                                +                                 ++e++dl+++ll++d+++Rlta+e++ hp++
  sp|P15790|CSK21_YEAST 295 KYLLKYEitlprefydmdqyirkpwhrfindgnkhlsgNDEIIDLIDNLLRYDHQERLTAKEAMGHPWF 363
                            5554444588999999****************************************************9 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (372 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.02
# Mc/sec: 4.91
//
Query:       sp|P21333|FLNA_HUMAN  [L=2647]
Description: Filamin-A OS=Homo sapiens GN=FLNA PE=1 SV=4
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------

   [No hits detected that satisfy reporting thresholds]

Domain annotation for each model (and alignments):

   [No targets detected that satisfy reporting thresholds]

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (2647 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         0  (0); expected 0.0 (0.02)
Passed bias filter:                        0  (0); expected 0.0 (0.02)
Passed Vit filter:                         0  (0); expected 0.0 (0.001)
Passed Fwd filter:                         0  (0); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               0  [number of targets reported over threshold]
# CPU time: 0.00u 0.00s 00:00:00.00 Elapsed: 00:00:00.00
# Mc/sec: inf
//

HMMER is a great tool and can also analyse multiple sequence alignments to build the hmms if that's what you need. I strongly urge you to reads its documentation which is both informative and a surpisingly fun read (it includes such phrases as "Their format is “proprietary”, which is an open source term of art that means both “I haven’t found time
to document them yet” and “I still might decide to change them arbitrarily without telling you”").


In-sample and across samples normalized expression I want to get the expression data that is in-sample normalized like FPKM and also across samples normalized as obtained using DESeq2 or else.

What I am currently doing is that I first normalize the data across samples (using DESeq) and from the resultant expression I calculate the FPKM. Does it make sense or am I missing something here?

This is how I am doing (not the exact code but the idea)

dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), ~ condition)
dds <- estimateSizeFactors(dds)
norm_data <- counts(dds, normalized=TRUE)

foreach sample
{
    foreach transcript
    {
        FPKM = (Normalized read count * 10^9) / (transcript length * total mapped normalized read count)
    }

}

Yes, this is a standard way of obtaining RPKM/FPKM/CPM values for plotting. Not that you do not need to use a for loop for any of the computations in R. You have a matrix of normalized values and things like transcript length are constant across samples (at least unless you're using something like salmon...although then you'd have TPMs to begin with).


Un-normalize DESeq2 counts I obtained a matrix of RNA-seq count data that has been normalized by DESeq2's median of ratio method.
I know that DESeq2 wants to take in un-normalized counts, but I do not have access to those data.
How do I best proceed here if I want to perform DEG analysis using DESeq2?
I know I can always start from .fastq files, but that would be so much extra work.
I don't think I can un-normalize the data because I don't know the scaling factors.
Any suggestions?

You can read in the normalized count table and don't normalize the data, but my advice here is not to do that. The rounding of the normalized matrix introduces some noise, but I think the larger issue is how sure are you that the table you are working with, is exactly a count table of normalized counts from DESeq2 ?
It really doesn't take much time to align and get count table again, rather than going back and forth wondering whether you did the normalization correctly.
I can show below with an example how it will look like, first we make an example dataset and obtain the normalized counts:
library(DESeq2)
set.seed(111)
sz = runif(6,min=0.5,max=1.5)
x = makeExampleDESeqDataSet(sizeFactors=sz,m=6)
x = estimateSizeFactors(x)
ncounts = counts(x,normalize=TRUE)
x = estimateDispersions(x)

Now suppose we feed this into a new DESeq2 object, but before that you need to do is check most of the size factors are close to 1:
estimateSizeFactorsForMatrix(round(ncounts))
 sample1  sample2  sample3  sample4  sample5  sample6
1.039179 1.045991 1.048276 1.045872 1.032420 1.052006

If you see something above that suggest a deviation from 1, something is wrong from the normalization and do not even think about using it.
Below you can use the normalized counts and you can see the difference is not huge:
dds = DESeqDataSetFromMatrix(round(ncounts),colData(x),~condition)
sizeFactors(dds) = 1
dds = estimateDispersions(dds)

Compare the new count table:
head(counts(dds,normalize=TRUE))
      sample1 sample2 sample3 sample4 sample5 sample6
gene1       0       2       0       1       0       0
gene2       9       8      16       6       1       2
gene3      23      40       8      38      11      41
gene4       0       0       0       6       8       0
gene5       3       2       0       5       1       2
gene6       0       1       3      19       4       9

And your input:
head(ncounts)
        sample1    sample2   sample3    sample4   sample5   sample6
gene1  0.000000  2.2658464  0.000000  0.9310803  0.000000  0.000000
gene2  8.693480  8.3081033 16.251453  5.5864817  1.121781  2.086619
gene3 23.472395 40.0299524  7.584012 38.1742919 11.217808 40.689076
gene4  0.000000  0.0000000  0.000000  5.5864817  7.852465  0.000000
gene5  3.477392  1.5105642  0.000000  4.6554014  1.121781  2.086619
gene6  0.000000  0.7552821  3.250291 18.6216058  4.487123  9.389787

We can check the dispersion, for lowly expressed genes you see some difference, but on the whole difference
with(elementMetadata(x),plot(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#FFC93C80"))
with(elementMetadata(dds),points(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#DBF6E980"))

Again, I would say I agree with @ATpoint comments that if you can, start from the fastq again.. don't work with data you have doubts about.


How to specify resources for cluster in snakemake I am trying to run snakemake version 3.5.4 on my cluster, where jobs are partially defined by a common template build according manual and also defined though rule-specific parameters. But for some reason cluster configuration is not expanding wildcards (and command does). When I run

snakemake data/genome.fa.gz.bwt -p --jobs 10 --cluster-config cluster.json --cluster "bsub -J "{rule}.{wildcards}" -q {cluster.queue} -n {threads} -M {resources.M} -R {cluster.one_host} -o {cluster.output} -e {cluster.error}"

-J "{rule}.{wildcards}" -q {cluster.queue} -n {threads} -M {resources.M} -R {cluster.one_host} expands as expected. But -o {cluster.output} -e {cluster.error} is replaced by literally logs/cluster/{rule}.{wildcards}.out and logs/cluster/{rule}.{wildcards}.err, i.e. there is no expansion in config file (the config file is attached bellow).

I am puzzled since exactly this lines are in the documentation. What am I doing wrong? Why snakemake is not expanding wildcards in the config file?

Here is my cluster.json

{
    "__default__" :
    {
        "queue"    : "bgee",
        "one_host" : "\"span[hosts=1]\"",
        "name"      : "{rule}.{wildcards}",
        "output"    : "logs/cluster/{rule}.{wildcards}.out",
        "error"     : "logs/cluster/{rule}.{wildcards}.err"
    },
}

Here is relevant part of Snakefile :

rule index_reference :
        threads: 1
        resources: M=20000000
        input : "data/{reference}/genome.fa.gz"
        output : "data/{reference}/genome.fa.gz.bwt"
        shell : "scripts/index_genome.sh {input}"

This feature was added in version 3.6.0. Look at release notes.

{wildcards} normally expands to a something that looks vaguely like a dictionary. In other words, something like key1=value1,key2=value2. Whether you really want to use that as a log file name is up to you.

Let's use a very contrived example to demonstrate this. Here is the Snakefile:

rule all:
    input: "input1.human.txt"

rule blah:
    output: "{sample}.{reference}.txt"
    shell: "touch {output}"

I'm using your cluster.json, so with that let's use a contrived cluster submission program:

snakemake --cluster-config cluster.json --cluster "echo 'rule.wildcards: {rule}.{wildcards} cluster: {cluster} -o {cluster.output} -e {cluster.error} -n {cluster.name}' >> log && " --jobs 1

If you run that and look in log, you can see that {rule} is expanded as expected, {wildcards} is expanded to sample=input1,reference=human. The output and error parts from cluster.json are expanded correctly for me. Below are the contents of log:

rule.wildcards: blah.sample=input1,reference=human cluster: one_host="span[hosts=1]",output=logs/cluster/blah.sample=input1,reference=human.out,name=blah.sample=input1,reference=human,queue=bgee,error=logs/cluster/blah.sample=input1,reference=human.err -o logs/cluster/blah.sample=input1,reference=human.out -e logs/cluster/blah.sample=input1,reference=human.err -n blah.sample=input1,reference=human

I'm not a big fan of using {wildcards} in paths unless you're doing a parameter sweep, but to each their own. Note that the above results are from snakemake-4.3.1.


How is the t-statistic value calculated in GEO2R or Bioconductor? I was trying to calculate t-statistics using Python's scipy and numpy as np-
scipy.stats.ttest_ind_from_stats(
    NP.mean(x[['GSM1224991', 'GSM1224992']]),
    NP.std(x[['GSM1224991', 'GSM1224992']]),
    len(['GSM1224991', 'GSM1224992']),
    NP.mean(x[['GSM1224993', 'GSM1224994']]),
    NP.std(x[['GSM1224993', 'GSM1224994']]),
    len(['GSM1224993', 'GSM1224994']),
)

For example in the dataset given below, GSM1224991 & GSM1224992 are in first group and GSM1224993 & GSM1224994 form second group.
ID_REF,GSM1224991,GSM1224992,GSM1224993,GSM1224994
148000,6.436150368369428,6.499787040655854,4.430816798843313,3.0819099697950434

Here is the calculation given in GEO2R analysis for the same reference but while I calculate t test it is 0.029630766122806237 instead of -0.33768 something as shown below.
ID  adj.P.Val   P.Value t   B   logFC
148000  1   0.001426    -9.44   -0.33768    -2.513

How to calculate t-stat here? Or considering single row is wrong approach?

That is not trivial. GEO2R uses the Bayesian linear model-based framework limma https://bioconductor.org/packages/release/bioc/html/limma.html for the analysis. I doubt you can easily and meaningfully implement the relevant code in Python without a major effort and notable background knowledge. The user guide https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf has background information plus links to the relevant papers. Don't try to analyze things by hand, use expert software. If it is not significant, well, then it is not significant. A 2vs2 comparison is strongly underpowered, therefore lack of DEGs are not unexpected, especially for microarrays.


Counting a specific consecutive character with its occurrence position and length I have a sequence file and want to count consecutive character "N", with the tandem's position of occurrence and its length. Say I have a file named mySequence.fasta like this:

>sequence-1
ATCGCTAGCATNNNNNNNNNNNNNNCTAGCATCATGCNNNNNNATACGCATCACANNNNNNNNNCgcatATCAC

and anticipated output should be like this:

Position 12 N 15
Position 38 N 6
Position 56 N 9

Kindly help me to solve this by awk or sed providing my file name mySequence.fasta

I couldn't come up with a one liner but have a two liner instead...

To get the start positions of the N repeat:

grep N mySeq.fasta | grep -b -o N | awk -F: '{print $1}' |awk '$1!=p+1{print $1}{p=$1}'
12
38
56

To get the length of the repeats:

grep N mySeq.fasta | awk -F '[^N]+' '{for (i=1; i<=NF; i++) if ($i != "") print length($i)}'
14
6
9

With thanks to stack overflow.


Declaring an array of strings and integers I have two csv file - one is a list of strings and other has both strings and floating values. If the strings in both the files match I need to copy the complete row from the second file to a temporary array. I have declared a matrix with n rows and m columns and saved is as dataframe temp <- as.data.frame(matrix(nrow = n, ncol = m)). I am using for loop to iterate over all rows and strcmp to compare the strings. But when the data is written to temp the strings are stored as integers. How should I declare this array? Is it possible to write this loop without declaring an array in R?

When you read in your two CSV files, you can explicitly declare the variable types for each column:

> a <- read.csv(..., colClasses=c("character"), ...)
> b <- read.csv(..., colClasses=c("character", "numeric", "numeric", ...), ...)

You could also investigate as.is or stringsAsFactors, although this can cause problems where there are duplicate strings. Using colClasses is more explicit.

Read the help documentation in ?read.csv for a more complete explanation of these operators.

You can filter with the %in% operator in match:

> strings_of_interest <- as.vector(a$V1)
> b.filtered <- b[b$V1 %in% strings_of_interest,]

If I am reading your question correctly, perhaps you can then take b.filtered and turn that into the matrix you want.

Avoid use of for loops where you can. They can slow R down considerably. A good read on this subject is in the third chapter of Patrick Burns' The R Inferno.


GSEA - Gene Set Enrichment Analysis - in R studio I have been trying to conduct a classic GSEA using the broad institute software and can not figure out how to compile it in R for the life of me. I am aware that there exists other packages in R that are almost the same, such as fgsea however, I would like to use this one first and then compare it with these other methods.

I have already generated a ranked list for my contrast of interest and save it as a .rnk file. Here is what I tried to do:

gseaCallH <-
  paste("java -Xmx2G -cp C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7",sep="")

testGSEA <- system(command = gseaCallH, intern = TRUE, ignore.stderr = TRUE)

When I run the system command I get the following error:
had status 1

At first I thought maybe it had something to do with JAVA not being properly installed however, I have installed it and still nothing. Maybe the installation was not done correctly? I have tried running it directly from cmd (I'm running windows 10) and got the following error:

Error occurred during initialization of VM
Could not reserve enough space for 2097152KB object heap

Can anyone advise on how to proceed?

I have answered my own question however, it is specific to windows users.. If you get the same error from the cmd as I did, you must manually put in environment variables related to Java. Do as the user Mohamed did in this SO post: https://stackoverflow.com/questions/9303889/error-occurred-during-initialization-of-vm-could-not-reserve-enough-space-for

Following this, you must avoid using paths such as C:/absolute_path and instead use /absolute_path, if you do the former, you will get the exception: java.net.MalformedURLException: unknown protocol: c

In the end, after manually entering the environment variables and running the following command: java -Xmx2G -cp /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file /Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx /Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out /Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7 it has worked.

Where the command can be executed directly from the command line or in R using the system function.

Note: A '#' separated by a space must be present at the beginning of your .rnk file for the program to work. Therefore, .rnk files must have the following format:

# Genes  Coeffecients
TSPAN6    0.2345676
etc. for every gene.


get gene lines from gtf file I would like to retrieve gene lines from a GTF file for which I only have exons & transcripts lines (output from Cufflinks) and alternative splicing possible. I need gene lines for compatibility with a pipeline dealing with GTF in Ensembl format.

Ideally, I would have gene lines representing the longest transcript - but I am open to discussion about that. The aim would be to separate the genome in genic & intergenic portions - boundaries being defined by start-end coordinates of the genes.

Sample input (e.g. gene "CUFF.105" with 5 transcripts):

cat subset_105.test | awk '$3=="transcript"{print $0}'

tig00000046 Cufflinks   transcript  26170   42766   202 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.3"; FPKM "0.2320304094"; frac "0.197788"; conf_lo "0.164885"; conf_hi "0.299176"; cov "5.188855";
tig00000046 Cufflinks   transcript  26170   42766   266 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.1"; FPKM "0.3061260755"; frac "0.249779"; conf_lo "0.228860"; conf_hi "0.383392"; cov "6.845843";
tig00000046 Cufflinks   transcript  26170   39469   470 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.2"; FPKM "0.5403628578"; frac "0.161731"; conf_lo "0.372512"; conf_hi "0.708214"; cov "12.084038";
tig00000046 Cufflinks   transcript  28928   39469   1000    +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.4"; FPKM "1.1485378983"; frac "0.233911"; conf_lo "0.853387"; conf_hi "1.443688"; cov "25.684547";
tig00000046 Cufflinks   transcript  29614   42766   181 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.5"; FPKM "0.2087076495"; frac "0.156792"; conf_lo "0.137112"; conf_hi "0.280304"; cov "4.667292";

What it came to my mind is:

sort the file by $4 (numeric: n) and $5 (reverse: nr) (sort -k4,4n -k5,5nr)
take the 1st line
remove unnecessary fields from $9

But this would fail for transcripts annotated in the reverse strand, and I won't be able to deal (how to fill) the $6  (score - a floating point value).

Any ideas of an existing method that could help me retrieve those gene lines?

Options to retrieve those lines directly from Cufflinks would also be accepted (I didn't find a suitable parameter in the manual).

Assuming none of the genes have duplicate entries across chromosomes (or strand or other biologically implausible things) in R one can:

library(rtracklayer)
library(GenomicRanges)
gtf = import.gff("file.gtf")
grl = split(gtf, gtf$gene_id)  # This produces a GRangesList
grl = endoapply(grl, function(x) {
    foo = x[1]
    foo$type = 'gene'
    start(foo) = min(start(x))
    end(foo) = max(end(x))
    return(c(foo, x))
    })
gr = unlist(gr)

You can then either write gr out to a file (either with export() or just printing the rows).


How to allow "half windows" in biopython? I am using Biopython to examine protein sequences including the Kyte and Doolittle hydrophobicity. However, for short sequences, this returns an empty value; the value is not allowed if the window does not have 20 values.

from Bio.SeqUtils.ProtParam import ProteinAnalysis

window_length=20
edge=1
kyte = {'A': 1.8, 'R': -4.5, 'N': -3.5, 'D': -3.5, 'C': 2.5, 'Q': -3.5, 'E': -3.5, 'G': -0.4, 'H': -3.2, 'I': 4.5,
        'L': 3.8, 'K': -3.9, 'M': 1.9, 'F': 2.8, 'P': -1.6, 'S': -0.8, 'T': -0.7, 'W': -0.9, 'Y': -1.3, 'V': 4.2}

full_sequence_analysis = ProteinAnalysis("CYGTGAVYLLGLGIGGFS")
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)

Is there are a built in Biopython option in this for "allow half windows/short windows" for incomplete windows? There is nothing I can see to suggest so in the docs.

You can choose between a window_size of 20 for sequences longer than 20 AAs, and the length of the sequence for those shorter::

seq = "CYGTGAVYLLGLGIGGFS"
window_length = min(len(seq), 20)

full_sequence_analysis = ProteinAnalysis(seq)
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)


How do you download stuff from NCBI fast? I have 2 .csv files that contains a list of accession codes. For example for this experiment a .csv file will contain:

SAMN10870244
SAMN10870243
SAMN10870242
SAMN10870241
SAMN10870240
SAMN10870239
SAMN10870238
SAMN10870237
SAMN10870236
SAMN10870235
SAMN10870234
SAMN10870233

Then I run the following

for i in $(cat *.csv); do echo ${i}; /usr/bin/sratoolkit/bin/fastq-dump --gzip ${i} & done

Some of my downloads takes a lot of time days. Is there a faster way to get stuff from NCBI?
Are there mirrors of NCBI to accelerate the process?
Why is it so slow?

Answers to this question can be found on biostars. To give you a short summary:

Faster way might be to use the parallel-fastq-dump, as suggested in this answer. I never tested that tool though, my own experience is that prefetch is more stable than fastq-dump command.
In that same post it seems that the Japanese mirror seems faster, don't know if that is still valid.
The reason this is never optimized might be that downloading from SRA is not a task you'll have to do repeatedly, as described here in this post.


fixing chromosome with new contig I want to fix a chromosome of around 50-100Mb (chr.fasta) for a region that is missing a few kilobases of data. I have the missing sequence with a window upstream and downstream that matches perfectly to the reference (contig.fasta).

What's the best way to produce a new chromosome, fixed.fasta that contains the new data from contig.fasta?

See depiction below:

chr.fasta    |---------------------------------------------------   -------------|
contig.fasta                                                   |-===-|

Final intended result:

fixed.fasta  |---------------------------------------------------===-------------|

Any ideas?

I have not received any answer to the question, so I created a small tool to fix the file myself.

I am assuming as per the question that there are perfect matches upstream and downstream to the new sequence, and that anything in between can be deleted/updated by the new sequence.

Given the first_n and last_n bases of the diagram in the question, and the entire patch sequence, we can create a single-line string of the reference with seqtk seq -l9999999 - and then edit it with the line below (in Perl regex):

$chromosome_str =~ s/$first_n.+$last_n/$patch/i

This creates the new string that can be used as the new fasta file.


Split fasta file based on groups in header information and output as separate files I have a fasta file containing the sequence of a gene across different species. In total there are around 900 samples and 12 species. (Each sequences is over multiple lines and longer than 100bp.)
My fasta file looks like:
>Species-1-samplenameA
CTATCCTTAAACGCATATCTCGCACAGTAACTCCCCAATATGTGAGCATCTGATGTTGCCCGGGCCGAGTTAGTCTTGTGCTCACGGAACTTATTGTATG
>Species-2-samplenameB
AGTAGTGATTTGAAAGAGTTGTCAGTTAGCTCGTTCAGGTAATGGTTCCTCACACTACGTCAAAATAAGAGAGCGGTCGTGACATTATCCGTGATTTTCT
>Species-3-samplenameC
CACTACTATCAGTACTCACGACTCGATTCTGCCGCAGCCACGTATCGCCAGAAAGCCAGTCAGCATTAAGGAGTGCTCTGGGCAGGACAACTCGCATAGT
>Species-3-samplenameD
GAGAGTTACATGTTCGTTGGGCTCTTCCGACACGAACCTCAGTTGGCCTACATCCTACCTGAGGTCTGTGCCCCGGTGGTGAGAAGTGCGCATTTCGTTC

I want to split this file in one fasta file per species.
I think it’s possible to use the awk function for this but I’m stuck. Does anyone have a script/code that might help me?
Thanks a lot.

A simple Biopython solution- iterate over the sequence records, identify the species, open the file handle using append mode to ensure no data is overwritten, and write the record:
from Bio import SeqIO

for record in SeqIO.parse("myfile.fa", "fasta"):
    species = record.id.split('_')[0]
    with open(f"{species}.fa", "a") as f:
        SeqIO.write(record, f, "fasta")


Are there databases to annotate non-coding mutations? I would like to ask if there are some known databases that allows to annotate non-coding mutations in a given region for the human genome?

Preferably, if the database is downloadable so I am able to automatize the processes of annotation.

VarSome allows you to annotate any variant, including non coding ones. While ACMG classification will not be available if the variant isn't in a gene's transcript (the ACMG criteria are only applicable to variants in transcribed regions), non-coding variants in transcripts will be fully annotated. For example, if you look up variant chr17-41197559-G-T which is a non-coding variant found in the UTR of one of the transcripts of BRCA1, you will get a lot of information (probably best to click on the link, the image has a lot of data):

For non-coding variants that don't fall in any transcript, you can still get frequency data, when available, as well as DANN pathogenicity predictions, GERP score, an overview of the genomic regions and other variants and/or srtuctural variants that can be found around the region, and any community annotations that may have been submitted by other users.

There is also an API you can use to automate the annotation, and VarSome Pro that lets you annotate an entire VCF file, but those are both paid services.

Disclaimer: I work for the company behind VarSome.


What's the meaning of having different alignments with the same qname in a bam with no duplicates? I have removed duplicates from my bam file (not secondary, nor suplementary alignments), and then I have looked for some repeated qnames (to assess the failure of a program) and I've found a lot of them. All of them has TLEN equal to 0.

As an example: this is part of a bed file that a generate from the bam.

CHR  Start   End       QNAME                                 PNEXT  TLEN
12 86945793 86945944 A05325:240:SS348DSZX:2:5105:12246:3790 90044835 0
4 168461445 168461485 A05325:240:SS348DSZX:2:5105:12246:3790 90044835 0

And those are these reads in the original bam:

A05325:240:SS348DSZX:2:5105:12246:3790  2177    4       168461445       0       111H40M 11      90044835        0       ATAAAATTAAAAAAATAAAAAAATTAAAAAAAAAAAAAAA        :,F,F,,FF:FF,F,,F,FFFF,,F,F,:,:,FFFFFFFF     NM:i:2  MD:Z:28T8G2     MC:Z:151M       AS:i:32 XS:i:31 RG:Z:HM327DSXX.3        SA:Z:12,86945793,+,20S37M94S,1,1;
A05325:240:SS348DSZX:2:5105:12246:3790  129     12      86945793        1       20S37M94S       11      90044835        0       ATGGTGAAAAAACAAATTATTTAAAAAAATCAAAAAAAAAAAAAAAACCAAAAAACATAAAAAAAAGGGCAGTAATTTTATTCATGGGGTCAAACAAAAAAAAAGCATAATATAAAATTAAAAAAATAAAAAAATTAAAAAAAAAAAAAAA      FF,,,:,F,F,,,,F,,,,,,,F,FFFF:,:FFFF:FFF,FF::F,:,,FFFF,,,F,FF::,,F,::,,,,F,,F:FF,,,,,,,F::F,F,,,,F,FFFFF:,F,,F,,:,F,F,,FF:FF,F,,F,FFFF,,F,F,:,:,FFFFFFFF      NM:i:1  MD:Z:9G27       MC:Z:151M       AS:i:32 XS:i:29 RG:Z:HM327DSXX.3        SA:Z:4,168461445,+,111S40M,0,2;

This are the commands used to remove duplicates:

samtools view -F 1024 -bo filtered.bam original.bam

The commands used to generate this bed were:

samtools view sample.bam | awk '{ print $3, $4, $4+length($10), $1, $8, $9}' > sample.bed

Why do I present those alignments?

I think that you are dealing with two different phenomena:

1) TLEN == 0 reads

From the BAM spec:

  TLEN: signed observed Template LENgth. If all segments are mapped to the same reference sequence, the absolute value of TLEN equals the distance between the mapped end of the template and the mapped start of the template, inclusively (i.e., end − start + 1).14 Note that mapped base is defined to be one that aligns to the reference as described by CIGAR, hence excludes soft-clipped bases. The TLEN field is positive for the leftmost segment of the template, negative for the rightmost, and the sign for any middle segment is undefined. If segments cover the same coordinates then the choice of which is leftmost and rightmost is arbitrary, but the two ends must still have differing signs. It is set as 0 for a single-segment template or when the information is unavailable (e.g., when the first or last segment of a multi-segment template is unmapped or when the two are mapped to different reference sequences).

Note that in each case the mate read appears to be mapped to a different chromosome/contig. So, the mapper can't estimate insert size in these cases, so it just says zero.

2. Multiple reads with same qname

There are various reasons why the same qname can appear multiple times. Again from the BAM spec, one such explanation is that it is a chimeric alignment:

  Chimeric alignment: An alignment of a read that cannot be represented as a linear alignment. A chimeric alignment is represented as a set of linear alignments that do not have large overlaps. Typically, one of the linear alignments in a chimeric alignment is considered the “representative” alignment, and the others are called “supplementary” and are distinguished by the supplementary alignment flag. All the SAM records in a chimeric alignment have the same QNAME and the same values for 0x40 and 0x80 flags (see Section 1.4). The decision regarding which linear alignment is representative is arbitrary.

You may have to do some digging and looking at the different flags to decide what exactly is going on, but my bet is that they are supplementary alignments.


Question: Bismark methylation for non-uniquely mapping BS-seq reads I am looking for the DNA methylation, specifically inside small RNAs; piRNA. I am mapping BS-seq reads on the genome with bismark. From the user guide of bismark, it maps reads uniquely on the genome, which mean it reports first read with best score and discard rest of the alignment and if a read mapping to multiple location with same score it will discard all such reads.
Typical command for obtaining uniquely mapped reads from bismark
~/Software/Bismark-0.22.3/bismark GenomePATH/Original_Genome/  --bowtie2 -p 10 -1 Methylation_reads/Sample1_R1.fq.gz -2 Methylation_reads/Sample1_R2.fq.gz --bam -o Sample1

Since, I am interested in piRNA only which may have multiple source in the genome (one piRNA coded from different location of genome), I wanted to take all mapped locations in consideration. Can you suggest what changes I have to make in bismark command. After reading bismark user guide I am confused because I could not able to find relevant parameter.
Also, please let me know if this approach is correct, if not what may be best approach for this.

I'm not certain about your piRNA approach, maybe others will be able to help with that question.
As for multi-mapped reads, from Bismark's documentation you can find the --ambig-bam flag:

For reads that have multiple alignments a random alignment is written out to a special file ending in .ambiguous.bam. The alignments are in Bowtie2 format and do not any contain Bismark specific entries such as the methylation call etc. These ambiguous BAM files are intended to be used as coverage estimators for variant callers.

It doesn't look like Bismark will return all of the multiple alignments, just a randomly selected one for a given read or read pair (this probably results from the way Bowtie2 handles multiple alignments within Bismark).


What is deep sequencing? People talk about deep sequencing. Is there any way to calculate how deep the sequencing is ? What should be the optimum depth to get reliable data ?

I am doing whole genome sequencing of a virus genome which is 10 kb long. I got 80000 reads from Illumina sequencing. Is there any way to tell about my sequencing depth ?

I found a post useful for this topic.

It explains the difference of coverage and depth.  It also has a useful explanation on how to calculate coverage and depth.

Here is a copy of what the link says just encase the post is removed:

  Depth of coverage

  How strong is a genome "covered" by sequenced fragments (short reads)?

  Per-base coverage is the average number of times a base of a genome is
  sequenced. The coverage depth of a genome is calculated as the number
  of bases of all short reads that match a genome divided by the length
  of this genome. It is often expressed as 1X, 2X, 3X,... (1, 2, or, 3
  times coverage).

  Breadth of coverage

  How much of a genome is "covered" by short reads? Are there regions that are not covered, even not by a single read?

  Breadth of coverage is the percentage of bases of a reference genome
  that are covered with a certain depth. For example: 90% of a genome is
  covered at 1X depth; and still 70% is covered at 5X depth.

SAMtools: get breadth of coverage

Get breadth of coverage of a reference genome.

How many bases of my reference genome are covered by my short reads sample at a level of 5X or higher?

 # set shell variables to your file-names

 REF_GENOME_FILE=my_ref_genome.fna
 WGS_SAMPLE=my_short_read_sample.fastq
 MIN_COVERAGE_DEPTH=5

Map short reads against reference genome (bowtie2)

# create bowtie2 index database (database name: refgenome)

bowtie2-build ${REF_GENOME_FILE} refgenome
ls
refgenome.1.bt2
refgenome.2.bt2
refgenome.3.bt2
refgenome.4.bt2
refgenome.rev.1.bt2
refgenome.rev.2.bt2

# map reads and sort bam file
bowtie2 -x refgenome --no-unal -U ${WGS_SAMPLE} -S - -p 12 | \
samtools view -bS - | \
samtools sort -m 5G - mapping_result_sorted.bam

-U ${WGS_SAMPLE} consider all reads as unpaired reads
-S -             bowtie output in SAM format, write to standard out
-x refgenome     use reference genome bowtie2-database
--no-unal        suppress SAM records for reads that failed to align
-p 12            use 12 processors

Get breadth of coverage (SAMtools)

# create samtools index
samtools index mapping_result_sorted.bam

# get total number of bases covered at MIN_COVERAGE_DEPTH or higher
samtools mpileup mapping_result_sorted.bam | awk -v
X="${MIN_COVERAGE_DEPTH}" '$4>=X' | wc -l
32876

# get length of reference genome
bowtie2-inspect -s refgenome | awk '{ FS = "\t" } ; BEGIN{L=0}; {L=L+$3};
END{print L}'
45678

Result: breadth of reference genome coverage
total number of covered bases: 32876  (with >= 5X coverage depth)
Depth of coverage (average per-base coverage): 0.719 X (32876 ÷ 45678) (total number of covered bases divided by reference genome length)
percent: 71.9%  (0.719 × 100) (coverage breadth in percent)


How to print alternative allele from VCF with htslib c++ I would like to print the alternative allele from a vcf file which i load in with htslib.
The description of some functions:
https://github.com/samtools/htslib/blob/develop/htslib/vcf.h
I have used bcf1_t objects by creating a header and reading the information.
htsFile *test_vcf = NULL;

  // creates header
  bcf_hdr_t *test_header = NULL;

  // initialize and allocate bcf1_t object
  bcf1_t *test_record = bcf_init();

  test_vcf = vcf_open("Testfile.vcf", "r");

  // returning a bcf_hdr_t struct
  test_header = bcf_hdr_read(test_vcf);

  while(bcf_read(test_vcf, test_header, test_record) == 0){
    std::cout << "chr" << bcf_hdr_id2name(test_header, test_record->rid) << ":" <<test_record->pos+1 << std::endl;
  }
  bcf_hdr_destroy(test_header);
  bcf_destroy(test_record);
  bcf_close(test_vcf);

But im not able to extract the reference or alternative allele information using this structure bcf1_t.
As far as i can read its perhaps possible with the structure and function:
bcf_dec_t
bcf_unpack()

But i am not able to get it to work since all the other functions seems limited to the bcf1_t struct.
Can anyone help with extracting the REF and ALT sequence?
Thanks for your time

see the code in https://github.com/samtools/htslib/blob/develop/vcf.c#L3104
 bcf_unpack((bcf1_t*)v, BCF_UN_ALL);// if you need genotypes too, might use another constant to avoid decoding all the structure
 (...)
 kputc('\t', s); // ALT
 if (v->n_allele > 1) {
        for (i = 1; i < v->n_allele; ++i) {
            if (i > 1) kputc(',', s);
            kputs(v->d.allele[i], s);
        }
    }


Assembly by stringtie I run this cmd

./stringtie G1_sorted.bam -B -o G1.gtf -G Triticum_aestivum.IWGSC.42.gtf -p 4  -C G1.refs.gtf -A G1.abund.tab

Error is:

-WARNING: no reference transcripts were found for the genomic sequences where reads were mapped!
Please make sure the -G annotation file uses the same naming convention for the genome sequences.

How can I resolve it?

Update:

I have downloaded both (GTF and Refernce Genome)files from Ensemble. Can you confirm me that whether theu have same chromosomes naming? So that i can perform my mapping  task of GTF File with Reference Genome which usually i remain fail to do due to difference in chromosomes naming.                                                           Here The Reference genome file is shown as:

 >1A dna:chromosome chromosome:IWGSC:1A:1:594102056:1 REF
CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCCTAACCCTAAACCCTAACCAAAAC

Where as The annotation file (gtf file ) is as

    #!genome-build IWGSC
#!genome-version IWGSC
#!genome-date 2018-07
#!genome-build-accession GCA_900519105.1
3B      IWGSC   gene    212892  214491  .       -       .       gene_id "TraesCS3B02G000100"; gene_source "IWGSC"; gene_biotype "protein_coding";
3B      IWGSC   transcript      212892  214491  .       -       .       gene_id "TraesCS3B02G000100"; transcript_id "TraesCS3B02G000100.1"; gene_source "IWGSC"; gene_biotype "protein_coding"; transcript_source "IWGSC"; transcript_biotype "protein_coding";

Can it be confirmend that these both files are same ? If not same with chromosome naming convention , then please suggest me how can I change it ?

This seems to be a common problem. As @llrs pointed out, Ensembl, UCSC(genome browser and associated data) and NCBI records are not directly interchangeable. The ambiguity refers to annotation as well as the reference genome. The reference genome is identical, except for the chromosome naming convention, meaning a reference from Ensembl might not contain CHR prefix, whereas it will be present in one downloaded from UCSC.

To prevent such conflicts in downstream tools, it is good to stick to the reference and annotation from the same repository throughout the analysis and map the genes to gene names or whatever towards the end of the investigation.

You can download all the associated data for T. aestivum which I beleive is the organism you want to study from Ensembl plants FTP download.


STAR quantMode vs RSEM vs Kallisto I recently discovered this Snakemake pipeline for RNASeq that uses STAR's quantMode to quantify gene expression for DESeq2 differential expression analysis.

In the past I've always seen the workflow as:

STAR alignment
RSEM/HTSeqCount quantification
DESeq2 analysis

More recently I've also seen:

Kallisto/Salmon pseudo-alignment quantification
DESeq2/Sleuth analysis

I tried searching online for differences between STAR quantMode vs other quantification algorithms but couldn't find many details. What are the benefits and drawbacks of using STAR quantMode vs RSEM/Kallisto/Salmon?

STAR quantMode (GeneCounts) essentially provides the same output as HTSeq-Count would, ie. number of reads that cover a given gene. This is the most simple measure of expression you could get from RNA-seq data. Kallisto and Salmon utilize pseudo-alignment to determine expression measures of transcripts (as opposed to genes). RSEM uses some algorithm to determine isoform fractions, but still based on alignments. (Note that I don't have any personal experience using RSEM.)

It ultimately depends on the research question which method will be most appropriate for you.

read counts

pros:

relatively easy to understand concept-wise
can be used to generate many different expression measures (FPKMs, TPM, CPM, etc.) and is, thus, more versatile in down-stream analyses

cons:

you can't distinguish isoforms (from the counts themselves)
requires alignment and, thus, more time to obtain

salmon/kallisto:

pros:

faster (due to lack of alignment step)
quantifies (known) isoforms

cons:

expression values are more abstract, less easy to explain
gene-level analyses will require some extra work

RSEM

Likely lies somewhere in between, still requiring alignments will make it slower. But it does give you expression values for isoforms.


Multi Omics comparison RNA seq : Condition specific gene signature The authors of an article are comparing different cell type and the organ type and they show the ontology enrichment of the RNA-seq:

Cell-type-specific and organ-specific marker genes were identified
using a two-step procedure. First, we performed pairwise differential
expression analysis for each of the three type of structural cells,
comparing them across organs. Second, based on the resulting pairwise
comparisons, we identified for each organ those genes that were
upregulated in each cell type compared to at least five other organs.
Pairwise comparisons between organs were performed using the limma
package in R, separately for each of the three structural cell types.
Significantly differential genes were selected based on statistical
significance (adjusted P < 0.05), average expression (log2(CPM) > −1)
and sequencing coverage (median number of reads greater than 10 in the
group with stronger signal). On the basis of these pairwise
comparisons, we counted the total number of times each gene was
upregulated in a specific organ compared to all other organs. Genes
that were upregulated in comparison to five or more other organs were
selected as marker genes of the corresponding organ.

Now how is this part done ? "Genes that were upregulated in comparison to five or more other organs were selected as marker genes of the corresponding organ.*"
The figure for this results is this:

From their methods section:

On the basis of these lists of marker genes and receptor–ligand pairs, we inferred potential cell–cell interactions for all pairs of one structural cell type and one haematopoietic immune cell type, quantifying the enrichment for known receptor–ligand pairs among all pairs of marker genes between the structural cell type and the haematopoietic cell type (Fig. 2a). First, we counted all pairs of marker genes for each pair of structural and immune cell types. Second, we calculated the fraction of these gene pairs that were annotated as receptor–ligand pairs. Third, we tested whether this fraction was greater than the fraction of annotated receptor–ligand pairs across all pairs of genes. Fisher’s exact test was used to obtain P values and odds ratios, as implemented by the function ‘fisher.test’ in R. P values were adjusted for multiple testing using the Benjamini–Hochberg method, as implemented the function ‘p.adjust’ in R. Finally, significantly enriched pairs (adjusted P < 0.05) of structural and immune cell types were connected by edges to generate a graph of cell–cell interactions.

In addition to testing for significant enrichment (as described in the previous paragraph), we calculated aggregated receptor and ligand gene signatures (Fig. 2c) based on a manually curated list of immune-related receptors and ligands (Supplementary Table 4). These gene signatures were derived by normalizing expression values for each gene and averaging across all genes in a given set of biologically related receptors and ligands.

The actual values plotted then seem to be from EnrichR, since they mention that elsewhere in the methods.


Any way to quantify the variation of genes that expressed in Affymetrix expression data? am experimenting preprocessed Affymetrix microarrays expression data matrix (Affymetrix probe-sets in rows (32830 probesets), and RNA samples in columns (735 samples)) for my downstream analysis. Here is how the data looks like:

> dim(eset_HTA20)
[1] 32830   735
## load assayData
HTA20_rma <- load("HTA20_RMA.RData")

> eset_HTA20[1:10,1:3]
             Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04
1_at                6.062215        6.125023        5.875502
10_at               3.796484        3.805305        3.450245
100_at              5.849338        6.191562        6.550525
1000_at             3.567779        3.452524        3.316134
10000_at            6.166815        5.678373        6.185059
100009613_at        4.443027        4.773199        4.393488
100009676_at        5.836522        6.143398        5.898364
10001_at            6.330018        5.601745        6.137984
10002_at            4.922339        4.711765        4.628124
10003_at            2.689344        2.771010        2.556756

objective:

I am interested in to reduce the dimension of the dataset by tossing off low expressed in genes in the experiment. To do so, I need to quantify the variation of the gene that expressed. I checked the possible approach and using the coefficient of variation (cv) could be one of the approaches I could try. It is not intuitive for me how to quantify the variation of genes that expressed.

Can anyone point me out how to lay out this procedure on the above expression data? Is there any Bioconductor package to do this kind of task? any idea to make this happen? thanks

You can use the following code to calculate the coefficient of variation:

# expr is your expression matrix.
SD <- apply(expr, 1, sd)
CV <- sqrt(exp(SD^2) - 1)

It might be implemented in some package but it is so brief that you can write it again yourself. Then you can filter out those that are below certain percentage of the distribution of CV (like 0.1)

expr[CV > quantile(CV, probs = 0.1), ]


What to use to edit RNA alignments? I have many alignments from Rfam Database, and I would like to edit them.
I saw that many tools are used for Protein sequence alignments, but there is something specific to edit RNA alignments ?

e.g. Stockholm Alignment of Pistol (just few entries).

FP929053.1/1669026-1668956   AGUGGUCACAGCCACUAUAAACA-GGGCUU-UAAGCUGUG-AGCGUUGACCGUC----------ACAA-----CGGCGGUCAGGUAGUC
AFOX01000025.1/1981-1912     ACUCGUCUGAGCGAGUAUAAACA-GGUCAU-UAAGCUCAG-AGCGUUCACCGGG----------AUCA------UUCGGUGAGGUUGGC
HE577054.1/3246821-3246752   ACUCGUCUGAGCGAGUAUAAACA-GGUCAU-UAAGCUCAG-AGCGUUCACCGGG----------AUCA------UGCGGUGAGGUUGGC
CP000154.1/3364237-3364168   GUUCGUCUGAGCGAACGCAAACA-GGCCAU-UAAGCUCAG-AGCGUUCACUGGA----------UUCG------UCCAGUGAGAUUGGC`
`#=GC SS_cons                 <<<<__AAAAA_>>>>-------..<<<<-.----aaaaa.----<<<<<<<<<..........____....._>>>>>>>>>-->>>>`
`#=GC RF                      acUCGUCuggGCGAguAUAAAuA..cgCaU.UAgGCccaG.AGCGUcccggcgg..........uUau.....uccgccgggGGUuGcg
//

I would suggest use RALEE—RNALignment Editor in Emacs. It can get for you the consensus secondary structure, you can move left/right sequences and their secondary structures (you can't do it in JalView!), and more.

It's an Emacs mode, so could be a bit hard to start off, but just try, you don't have to use all Emacs features to edit your alignments!

  The RALEE (RNA ALignment Editor in Emacs) tool provides a simple
  environment for RNA multiple sequence alignment editing, including
  structure-specific colour schemes, utilizing helper applications for
  structure prediction and many more conventional editing functions.

Sam Griffiths-Jones Bioinformatics (2005) 21 (2): 257-259.

Fig. You can move left/right sequences and their secondary structures (you can't do it in JalView!)


Software to produce a table of post-translational modifications from a peptide list Does anyone know if there is a program/library/script in R or Python that takes as input a list of proteins/peptides and a list of post-translational modifications (PTMs; like oxidation of methionine and acetylation of cysteine), and as output returns a matrix/list with the imported peptides and all possible locations of the selected modifications?

I don't know of a tool that can do exactly what you want. However, there are various tools that can predict specific classes of post translational modifications. For example:

Phosphorylation: GPS

Computational prediction of phosphorylation sites with their cognate protein kinases (PKs) is greatly helpful for further experimental design. Although ~10 online predictors were developed, the PK classification and control of false positive rate (FPR) were not well addressed. Here we adopted a well-established rule to classify PKs into a hierarchical structure with four levels.

Mitochondrial import: Mitoprot

MitoProt calculates the N-terminal protein region that can support a Mitochondrial Targeting Sequence and the cleavage site. A complete description of the method to make the prediction is available in: M.G. Claros, P. Vincens. Computational method to predict mitochondrially imported proteins and their targeting sequences. Eur. J. Biochem. 241, 779-786 (1996).

Acetylation: "PAIL: Prediction of Acetylation on Internal Lysines"

In this work, we present a novel online predictor for protein acetylation sites prediction of PAIL, Prediction of Acetylation on Internal Lysines. We have manually mined scientific literature to collect 249 experimentally verified acetylation sites of 92 distinct proteins. Then the BDM (Bayesian Discriminant Method) algorithm has been employed.

You can use tools like the above and then combine the results to get the sort of matrix you are asking for.


Executing PyMOL from a Shell script I'm trying to execute pymol from a shell script (and it is not working). I'm not executing the script on  PyMOL but even if I do this, it doesn't work neither. However, if instead of a script the commands are manually written in PyMOL, everything goes well.

Here is the code:

#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
' > scriptt.pml
pymol scriptt.pml

One of the problems, you will see, is that PyMOL does not remove the water molecules. Another issue is that the program does not exit. Do you have any idea what's going on there?

Your shell script is fine, and works as expected. The reason it hangs is that you aren't telling PyMol to exit after running your script, so it's still running and waiting for the next command. Just add a quit at the end and it should work:

#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
select prot, 1OYT and not ligand,
save lig.pdb, ligand
save prot.pdb, prot
quit
' > scriptt.pml
pymol scriptt.pml
obabel prot.pdb -xr -O prot.pdbqt
obabel lig.pdb -O lig.pdbqt
echo 'babel conversion is finished'


Convert SRA to FastA I'm trying to get the FastA files for some accessions (like NC_001416.1). I did not managed to find an FTP server or direct link to these files (I want to get it from command line with wget, not from a web browser). But I found an "equivalent" file in ftp://ftp.ncbi.nlm.nih.gov/sra/refseq/NC_001416.1

This file is an SRA file that should be processed with the SRA-toolkit. In particular, fastq-dump seems to be the tools of choice.

I ran

fastq-dump --fasta NC_001416.1

and I got a NC_001416.1.fasta file. But the file contains the sequence I'm looking for in chunks of 5000 bases. The IDs in the output file are:

>NC_001416.1.1 length=5000
>NC_001416.1.2 length=5000
>NC_001416.1.3 length=5000
>NC_001416.1.4 length=5000
>NC_001416.1.5 length=5000
>NC_001416.1.6 length=5000
>NC_001416.1.7 length=5000
>NC_001416.1.8 length=5000
>NC_001416.1.9 length=5000
>NC_001416.1.10 length=3502

I need that sequence under a single ID:

>NC_001416.1 length=48502

I cannot find the proper parameters to supply to fastq-dump to achieve the desired result.

Does anyone know how could I get the FastA file I'm expecting? How to process properly the SRA file or where I could download the FastA file directly.

Finally, I found an alternative to the SRA translation: a link that works! For those of you interested in knowing how to download FastA files from NCBI using an accession number, try the following link:

https://www.ncbi.nlm.nih.gov/search/api/sequence/${accession}/?report=fasta

Using wget to download the accession used as example:

wget -O NC_001416.1.fasta "https://www.ncbi.nlm.nih.gov/search/api/sequence/NC_001416.1/?report=fasta"


Which reference to use for read mapping for popular model organisms What is the "best" assembly for the popular model organisms:

human (GRCh37 and GRCh38 are obvious, I'd pick whatever bwakit uses)
mouse (GRCm37/GRCm38, OK)

but what about non-human/mouse ones?

fruit fly
zebrafish
E. coli
any other idea?

I will answer one of the points – E.coli.

TL;DR Bacteria, and in particular E.coli, are highly variable and there is usually no single best assembly. Large scale WGS studies should come with multiple assemblies for individual monophyletic clusters.

Long answer: Whereas a single reference sequence can make sense for human or mice (to a certain extent), bacteria are not the case. The variability between them even within a single species can be enormous and various lineages can contain very different genes. This is the reason why biologists usually distinguish a core genome consisting of the genes shared among a vast majority (typically 95%) of individuals within a given phylogenetic clade, (i.e., of those genes which are essential for this clade – imagine a CPU and RAM in a computer), and an accessory genome consisting of the other, unnecessary, genes – imagine a joystick in our analogy). The size of the accessory genome can vary a lot for different bacteria. For instance, Escherichia coli is said to have an open accessory genome (i.e., much bigger than the core genome; imagine a Raspberry-PI and all the extension kits) whereas Chlamydia trachomatis has a closed accessory genome (i.e., small compared to the core genome; imagine an ATM).

For species having a small accessory genome, a single reference might be sufficient. However, if a core genome is large, we can use either pan-genomes, or multiple references. The word pan-genome tends to be used for different things: the computational biologists usually use it as a representation of all genomic content in a certain phylogenetic clade (which can be mathematically modeled using finite automata, hence as graphs) whereas some biologists define it as a concatenation of all (core and accessory) genes from the clade (without considering SNPs, etc.). Unfortunately, methods for graph references (especially for read mapping) are still not sufficiently developed, therefore, researchers still have to resort to reference sequence-based methods and use either concatenated genes or multiple reference sequences.

When particular bacterial species are studied using massive whole-genome sequencing of many isolates, researchers usually cluster the isolates to monophyletic sequence clusters and infer a special reference sequence for each of these clusters (see, e.g., this example for the pneumococcus). All of the inferred references should be then published together with the paper or in a separate data paper (see, e.g., the sequences from the previous example).

Disclaimer: I am not a biologist. If anything in this post is too simplified or incorrect, please, leave me a comment and I will try to fix it. I tried to make this text accessible for my younger self.


How do I export my biom object to a file in biom format? I'm using the python API for the biom-format package. I read in two files in biom format and merged them, like so:

In [1]: import biom

In [2]: a = biom.load_table('./table-a.biom')

In [3]: b = biom.load_table('./table-b.biom')

In [4]: a
Out[4]: 51116 x 46 <class 'biom.table.Table'> with 185077 nonzero entries (7% dense)

In [5]: b
Out[5]: 37479 x 26 <class 'biom.table.Table'> with 100124 nonzero entries (10% dense)

In [6]: c = a.merge(b)

In [7]: c
Out[7]: 67022 x 47 <class 'biom.table.Table'> with 254555 nonzero entries (8% dense)

How do I now export this biom object (c) as a biom file? There are a few methods available, eg, .to_tsv(), but I don't see a way to save this object in biom format.

You can use the functions to_hdf5 and to_json from the biom-format package:

to_hdf5

with biom.util.biom_open('output.biom', 'w') as f:
    merged_table.to_hdf5(f, "example")

to_json

with open('output.txt', 'w') as f:
    f.write(merged_table.to_json("example"))


Genome assembly from error-prone reads I understand how to assemble genome from error-free reads. I implemented like this:

  Construct directed overlap graph with reads as vertices and edges as
  maximum overlap between two vertices. Edges represent the length of
  overlapping maximum substring of one read suffix with other read
  prefix. Find the path that visits all vertices only once in overlap
  graph. With small graph we can use greedy approach : Starting from
  first pattern always choose egde of maximum weight and store shortest
  common superstring spelled by this path in resulting string. Resulting
  string spelled by path is shortest common superstring of all reads,
  that is assembled genome.

For example

AACTAG
  CTAGAT

Gives the AACTAGAT as a result

Now let's say that two error-prone reads have at most one error (mismatch):

AACTAG
  CTCGAT

As you can see second read has one substitution in third character from 'A' to 'C'.

So how assembled genome should look like for the error-prone reads? Do I take A or C to resulting string? Or I'm missing something?

Note: I used short reads for example because of simplicity. I work with much larger number of reads which are 100 nucleotides long and have at most one substitution of character. So that's a big overlap graph. I know I can use De Bruijn graph, but can I use some greedy approach (like in error-free example)?

For short reads, the typical and the most widely used solution is to correct away sequencing errors before assembly. You can correct errors with k-mer spectrum, a trie or multi-alignment. There are many papers on this topic.

Error correction alone won't fix all sequencing errors. Remaining errors may lead to bubbles and tips in the overlap/de Bruijn graph. You can identify them by inspecting the local topology of assembly graph.