Consensus on how to select SNPs for creating a phylogeny when the full dataset is too large I have got a quite large VCF file (about 20k samples with 160k SNPs after filtering for quality etc.) and I would like to get a phylogeny for it. However, the whole dataset's too large for my computational resources (unless I go for fasttree which would be my last resort). Therefore, I need to subset the SNPs and there's plenty of possible ways to do this (dump all rare SNPs with MAF below a threshold, only use SNPs in coding regions, only use missense SNPs; just to name a few).
Is there a consensus (or at least some comparative work or recent review; I found nothing) on what method works best and gives the tree that's closest to using the full dataset?

The subset that gives you the tree closest to the one you obtain with the full dataset does not necessarily is the tree nearest to reality... It could be that you get a more realistic result using only partial data for multiple reasons. Your best chance is in my opinion, as it usually is, making lots of tests and contrasting hypotheses.


genotyping or variant calling in R Are there any good options for calling variants or even just a pileup from R? R is not ideal for this, but I'd like to integrate with other functions.
I found deepSNV::bam2R which roughly does what I am looking for. Are there any alternatives?

Samtools/bcftools/tabix is available in R via Rsamtools.
SEQprocess also provides wrappers for different tools such as GATK and bowtie.
More utility packages are noted here.


maftools tcgaCompare plot: x-axis labels getting clipped I have a {maftools} plot:
whatever I try to save that completely I see the down part of plot is trimmed like

I have tried both save as pdf and save as image options which did not differ
Any suggestion please?
Plot comes by this code
tcgaCompare(maf = laml,cohortName = "Cohort",rm_zero = FALSE,decreasing = TRUE,logscale = TRUE,primarySite = TRUE,col = c("grey30", "red"),medianCol = "black")

One way is to increase the outer margins of plot.
par(oma=c(7,0,0,0))
tcgaCompare(maf = laml, primarySite = TRUE)


How to get the product of a CDS I need the name of the protein in /product="protein_name" using bash commands. Beware, there is a lot of whitespace between lines.

FEATURES             Location/Qualifiers
     source          1..1266
                     /organism="Sarcophilus harrisii"
                     /mol_type="mRNA"
                     /isolate="Salem"
                     /db_xref="taxon:9305"
                     /chromosome="1"
                     /map="unlocalized"
                     /sex="female"
                     /country="Australia: Taronga Zoo, Sydney"
     gene            1..1266
                     /gene="UROS"
                     /note="Derived by automated computational analysis using
                     gene prediction method: Gnomon. Supporting evidence
                     includes similarity to: 10 Proteins, and 100% coverage of
                     the annotated genomic feature by RNAseq alignments,
                     including 5 samples with support for all annotated
                     introns"
                     /db_xref="GeneID:100933602"
     CDS             60..908
                     /gene="UROS"
                     /codon_start=1
                     /product="uroporphyrinogen-III synthase isoform X1"
                     /protein_id="XP_003755055.2"
                     /db_xref="GeneID:100933602"
                     /translation="MRIPTSFSRGAWWKMKVLLLKDAKDDDSGLDPYLQELGSCGFEA
                     TLIPVLAFEFLSLPSFSEKLSHPEAYGGLVFTSPRAVEALQLSLERAGQLEVWKKSLK
                     EKWRLKPAYVVGQATASLVKKVGLAAEGEKCGNGEKLAEYICSRETPATLPLLFPCGT
                     LKGEVLPKTLKQKGILLESLTVYQKVQHPRMQDSLKSYFSTQGVPAGVTFFSPSGLTY
                     SLPHIQELSGASFDRIKFAAIGPTTARALASAGIPVSCTAEKPCPQDLAAGLQKALQP
                     SGPSSS"

Right now I have this:

while read p;
do cat "${p}.gb" | tr "\"" " "| egrep -v "/db_xref= taxon"| egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk '{print($2)}END{print("")}'
done <XM_codes.txt

Which results in this:

GeneID:100933602
UROS
uroporphyrinogen-111
XP_003755055.2

My output should be:

GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2

If you change the delimiter to something besides a space in the tr command (tr "\"" "\t") then you can use that same delimiter in the awk command (-F "\t").  In my example below I used a tab, also I filtered the taxon prior to the tr so that you don't have to worry about that staying in sync if you change the delimiter.

 egrep -v "/db_xref=\"taxon"| tr "\"" "\t" | egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk -F "\t" '{print($2)}END{print("")}';

results:

GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2


Number of reference sequences in a SAM file From a single record, I can get the reference sequence ID from the field rID, but how can I get the total number of different reference sequences stored in the whole SAM file?
It can be simple as just looping through all records, but I need some other solution.

There are multiple ways. Here are two:

samtools idxstats ‹bamfile›
samtools view -H ‹bamfile› | grep '^@SQ'

Both commands give you a list of the references in the BAM file. To get their number, simply append | wc -l. But note that idxstats also outputs a final entry, *, which doesn’t correspond to any reference but to all unmapped reads. So subtract that from the number.

However, beware that while this lists all references associated with a BAM file, not every reference necessarily has records associated with it. So there might be no reads mapping to a given reference. This is the third column in the idxstats output. So, to get only those references with at least one mapped read, we can filter for those:

samtools idxstats ‹bamfile› | awk '$4 != 0' | head -n-1 | wc -l

Oh, and these (and other) commands work on BAM files, not SAM. The idxstats command additionally needs the BAM index, since it’s actually reading information from that index file, not the BAM (this also means that it only works on sorted BAM files). Working on SAM would require looping through the whole file, which is quite inefficient. You should generally be working with BAM rather than SAM, so this is hopefully not a restriction.

To do the same on a SAM file, you need to iterate over the file:

awk '{print $3}' ‹samfile› | sort -u | wc -l

Or:

awk '{refs[$3] = 1} END {for (ref in refs) print ref}' ‹samfile› | wc -l


Why does a missing label 11plex TMT shows up at almost 50% intensity compared to other labels? We have a test plate that was done with 10 of the 11plex TMT labels. So, the intensities of the 11th label should be almost zero, plus the leakage from the 10th label due to 13C content. Although we observe a substantial drop in the intensities, we see ratios up to 50 percent compared to the highest reference intensity, which is the maximum of the intensities of the other labels, if we sum up the intensities for all proteins. For individual proteins the intensity sometimes is higher than for all "real" samples. Only talking about proteins that are identified "Confidential". Is that something you would expect? Is this due to the normalization?
I was using Reporter 0.7.17 from command line without specifying a reference sample. Pretty much default configuration and based on the default protein report.

Fig: Ratios (intensities) reported by the program.

The developers answered along these lines:
This might be a normalization artefact. Reporter is normalizing the ratios so that the median is 1. So, if the label is almost not present the noise gets increased. To check that you can turn of the normalization. This function is currently not documented. The developers are working on it though.

GitHub issue


How can I obtain the percentage gene expression per identity class in Seurat as further processible numbers (e.g. matrix)? I am analysing my single cell RNA seq data with the Seurat package. I want to know if there is a possibilty to obtain the percentage expression of a list of genes per identity class, as actual numbers (e.g. a matrix) which I can write out to say an excel file. For example I want to be able to say gene X is expressed in 30% of cell A, 10% of cell B etc.

I can generate something comparable with:

DotPlot(object = my_object, genes.plot = "my_gene")

However the results are only graphic and I wish to have further processible numbers.

Furthermore:

AverageExpression(object, genes.use = "my_gene")

Produces expression values which I cannot transform to percentages.

I will be very grateful on any hints.

This can be solved like this:

library(Seurat)

my_genes <- c("gene1", "gene2", "gene3")

exp <- FetchData(object, my_genes)

matrix <- as.matrix(colMeans(exp  > 0))*100


Parallel processing of scripts that use obitools I have an obitools script (de Barba et al. 2016) that I would like to run faster. How would you run it in parallel to cut down on time?

illuminapairedend -r rawdata_scandinavia_R2.fastq rawdata_scandinavia_R1.fastq | tee rawdata_scandinavia.fastq | obiannotate -S goodAli:'"Alignement" if score>40.00 else "Bad"' | obisplit -t goodAli -p rawdata_scandinavia.
touch rawdata_scandinavia.Bad.fastq
touch rawdata_scandinavia.Alignement.fastq
touch rawdata_scandinavia.fastq
ngsfilter -t rawdata_scandinavia.ngsfilter -u rawdata_scandinavia.unidentified.fastq rawdata_scandinavia.Alignement.fastq > rawdata_scandinavia.filtered.fastq
obisplit -p MICROSAT.PCR_ -t experiment rawdata_scandinavia.filtered.fastq

Can someone describe what tools and steps one needs to take to run this job in parallel, given that the functions are not designed to be run in parallel out of the box?

De Barba, M., Miquel, C., Lobréaux, S., Quenette, P. Y., Swenson, J. E., & Taberlet, P. (2016). High-throughput microsatellite genotyping in ecology: improved accuracy, efficiency, standardization and success with low-quantity and degraded DNA. Molecular Ecology Resources, 1–16. https://doi.org/10.1111/1755-0998.12594

Your only option here is to first split the fastq files into small chunks, process each and then merge them back together. Alternatively, find something else that runs faster.


Converting a VCF into a FASTA given a reference with Python, R I am interested in converting a VCF file into a FASTA file given a reference sequence with Python or R.

Samtools/BCFtools (Heng Li) provides a Perl script vcfutils.pl which does this, the function vcf2fq (lines 469-528)

This script has been modified by others to convert InDels as well, e.g. this by David Eccles

./vcf2fq.pl -f <input.fasta> <all-site.vcf> > <output.fastq>

Given that converting a VCF into a FASTA given some reference is a somewhat standard procedure in bioinformatics, are there standard R or Python implementations used by the community?

You can whip up something quite easily in Python using pyfaidx, which allows you to access a FASTA file by genomic coordinates. So you'd just have to pull the coordinates from the VCF to get the sequence, then output in FASTA format with whatever unique identifier for each sequence you want to use.

Edit: The FastaVariant class is new in pyfaidx, I hadn't seen it before. It makes this even easier, as you won't have to manually mess with the VCF to get the positions. Something like the following should work:

consensus = FastaVariant('yourgenome.fasta', 'variants.vcf.gz', het=True, hom=True)

out = open("variants.fasta", "w")

for chrom in consensus.keys:
    for var in consensus[chrom].variants_sites:
        record = consensus[chrom][var-1:var]
        print(record, file=out)

out.close()


merge the matrix from computematrix of deeptools I have several matrix from computematrix of deeptools. I need to merge two of them using "computeMatrixOperations cbind -m input1.mat.gz input2.mat.gz -o output.mat.gz"
but I am running to error "/computeMatrixOperations.py", line 521, in cbindMatrices
    hm.matrix.matrix = np.hstack((hm.matrix.matrix, np.empty(hm2.matrix.matrix.shape)))"

Is there any other way to merge the matrixes?

You're trying to cbind two matrices of different sizes, so computeMatrix is throwing an error since it's unclear which row goes with which. You need to ensure that matrices have exactly the same number of rows and that the Nth row in each match (i.e., they should be bound together with cbind).


BioJava - Protein to DNA Using Biojava 5, I am trying to convert a sequence of Amino Acids to a Nucleotide sequence. Any idea how to do this?

Note that unlike the reverse (nucleotide sequence to protein sequence), which is uniquely defined given a reading frame, the mapping from protein sequence to nucleotide sequence is not unique (it is not a function in the mathematical sense). The method you use to select which of the many (exponential in input sequence length) possible sequences to output depends critically on what the question is.
If the goal is to design a synthetic gene that will produce a desired protein, then you need to optimize it for expression in the organism into which the gene will be inserted, subject to the constraint of the encoded amino acid sequence.
If the goal is to search for nucleotide sequences in one or more genomes that could give rise to a certain polypeptide, there is a much better way than to do reverse translation at all. In particular, most putative protein-coding open reading frames in most sequenced genomes have already been run through forward translation code and the expected amino acid sequence encoded in them can be found in databases. So, this problem is reduced to a substring search in these databases.


How to separate a microbial strain of MTB resistant and susceptible drugwise? Im trying to classify my list into three type: XDR - A strain must be resistant to:

MDR + CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN + AMYCACIN/CAPREOMYCIN/KANAMYCIN-> means ISONIAZID + RIFAMPICIN + any of these or two or three or all (CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN) + any of these or two or three or all (AMYCACIN/CAPREOMYCIN/KANAMYCIN)

PRE-XDR - A strain must be resistant to:

MDR + CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN -> means ISONIAZID + RIFAMPICIN + any of these or two or three or all (CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN)

MDR - A strain must be resistant to:

ISONIAZID + RIFAMPICIN
and if it is resistant to any one/two/three of these or all of these in
addition to ISONIAZID + RIFAMPICIN than again it is MDR:-
PYRAZINAMIDE, ETHIONAMIDE, STREPTOMYCIN, ETHAMBUTOL

SUSCEPTIBLE - A strain must be Susceptible to all drugs

I tried the code as follows for MDR:

tail -n +2 drug.txt |\
awk '{a[$1] = a[$1]","$2};END{for(i in a)print i"\t"a[i]}'|\
grep -P "\t,ISONIAZID,RIFAMPICIN$"

It gave me MDR starins correctly but than i had to first remove all Susceptible strains and applied it to Resistant one only.

But when i applied the code for Pre-XDR:

tail -n +2 drug.txt |\
awk '{a[$1] = a[$1]","$2};END{for(i in a)print i"\t"a[i]}'|\
grep -P
"\t,ISONIAZID,RIFAMPICIN|CIPROFLOXACIN|MOXIFLOXACIN|OFLOXACIN|LEVOFLOXACIN"

It is giving me the result as below:

SRR671821       ,CAPREOMYCIN,ETHAMBUTOL,ETHIONAMIDE,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN,STREPTOMYCIN
SRR671863       ,CAPREOMYCIN,ETHAMBUTOL,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN
SRR671822       ,ETHAMBUTOL,ISONIAZID,OFLOXACIN,RIFAMPICIN,STREPTOMYCIN
SRR671864       ,CAPREOMYCIN,ETHAMBUTOL,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN

Which is not correct for Pre-XDR.

So I'm looking for a the drugs if match than look for resistant and susceptible if these criteria match with any of the mention type (MDR,XDR...) it put those under that in a new file. Can anyone help?
my file druglist.sxls

MDR eg from the file attached is:

ERR038736   ISONIAZID   RESISTANT
ERR038736   RIFAMPICIN  RESISTANT
ERR038737   ISONIAZID   RESISTANT
ERR038737   RIFAMPICIN  RESISTANT
SRR1270146  AMIKACIN    SUSCEPTIBLE
SRR1270146  CAPREOMYCIN SUSCEPTIBLE
SRR1270146  CYCLOSERINE SUSCEPTIBLE
SRR1270146  ETHAMBUTOL  RESISTANT
SRR1270146  ETHIONAMIDE SUSCEPTIBLE
SRR1270146  ISONIAZID   RESISTANT
SRR1270146  KANAMYCIN   SUSCEPTIBLE
SRR1270146  PARA-AMINOSALICYLIC ACID    SUSCEPTIBLE
SRR1270146  PYRAZINAMIDE    RESISTANT
SRR1270146  RIFAMPICIN  RESISTANT
SRR1270146  STREPTOMYCIN    RESISTANT

Since this involves combining multiple conditions, it is easier (and certainly more readable if you ever need to come back to it later) to write a script to do this. For example, in perl:

#!/usr/bin/perl

use strict;
use warnings;

## This will hold the information for each strain
my %strains;

## These are the functions that will test each strain
sub isMDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}){
    foreach my $drug (keys(%thisStrain)) {
      next if $drug eq "ISONIAZID";
      next if $drug eq "RIFAMPICIN";
      next if $drug eq "PYRAZINAMID";
      next if $drug eq "ETHAMBUTOL";
      next if $drug eq "STREPTOMYCIN";
      ## discard this strain if it is resistant to any other drug
      if ($thisStrain{$drug}) {
        return 0;
      }
    }

    return 1;
  }
  return 0;
}

sub isPreXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}) &&
  (
   $thisStrain{"CIPROFLOXACIN"}  ||
       $thisStrain{"MOXIFLOXACIN"} ||
   $thisStrain{"OFLOXACIN"}  ||
       $thisStrain{"LEVOFLOXACIN"}
      )
     ) {
    return 1;
  }
  return 0;
}

sub isXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  isPreXDR($_strain) &&
      (
       $thisStrain{"AMYCACIN"} ||
   $thisStrain{"CAPREOMYCIN"} ||
       $thisStrain{"KANAMYCIN"}
      )
     ){
    return 1;
  }
  return 0;
}

sub isSusceptible{
  my $_strain = shift;
  foreach my $drug (keys(%{$strains{$_strain}})) {
if ($strains{$_strain}{$drug}) {
      return 0
    }
  }
  return 1
}

## Iterate over each line in the input file.
## For each iteration, the special variable $_
# will hold the contents of the line.
while (<>) {
  ## Skip first line
  next if $. == 1;
  ## remove Windows-style line endings
  s/\r//;
  ## remove trailing newline
  chomp;
  ## split the line on commas
  my ($strain, $drug, $status) = split(/,/);
  if ($status eq "SUSCEPTIBLE"){
$status = 0
  }
  else {
$status = 1;
  }
  $strains{$strain}{$drug}=$status;
}

foreach my $strain (keys(%strains)) {
  print "$strain\t";
  if (isMDR($strain)) {
print "MDR\n";
  }
  elsif (isXDR($strain)) {
print "XDR\n"
  }
  elsif (isPreXDR($strain)) {
print "PRE-XDR\n"
  }
  elsif (isSusceptible($strain)) {
    print "SUSCEPTIBLE\n";
  }
  else {
    print "OTHER\n";
  }
}

If you save that script as foo.pl, you can run it on your file with:

perl foo.pl drug.csv

Selected output from running this on the drug.csv file you posted is:

$ foo.pl drug.csv | grep -E 'ERR038288|ERR038736|ERR038751|SRS553990|ERR133879'
ERR038751   OTHER
ERR038288   SUSCEPTIBLE
SRS553990   XDR
ERR038736   MDR
ERR133879   PRE-XDR

Various details about this question were added in chat: https://chat.stackexchange.com/transcript/message/47821646#47821646


How would we test for significant changes in expression when many of the normalised counts are clustered together? I've been handed RNA-seq data with a lot of associated covariates.
This data has been put through the DESeq2 package and as a result normalised the data.
One of the transcripts of interest still has non-normal distribution for two genotypes. This is because many of the normalised counts are now clustered around a value of 4.

This is making it tricky to perform statistical tests because in order to perform parametric tests (in this case I want to perform a ANCOVA or multiple linear regression to include the covariates) the data must be normally distributed.
I also don't want to use a non-parametric test because I can't include all the covariates.
Should I just perform parametric tests ignoring this violation or do I have to just perform non-parametric on all my data for a fair comparrison?
Or is there another alternative option I am missing?
Thanks in advance

TLDR- just do the analysis in DESeq2 (it can handle multiple covariates)
If I've interpreted your question correctly there seem to be a few confusions here.
First off, for a multiple linear regression aka ANCOVA, one of the assumptions are the residuals are normally distributed - however this is almost never an issue in practice. The distribution of the data is not so important, you make the assumption that the data has come from a population which is continuous and normally distributed that's why you've chosen to do an ANCOVA. But, RNAseq produces count data - not continuous and bounded at zero. You need to be looking at using a Poisson GLM, or negative binomial GLM as done by DESeq2 (lookup overdispersion)
The normalisation done by DESeq is to account for different library sizes,it doesn't transform the data to make it conform to a normal distribution. You'd need to do some form of data transformation to do that, but don't it's unnecessary.
I'd suggest reading the DESeq vignettes and following their standard pipeline.
Hope that is of some help.


How to check whether all BAM read contain defined read groups? I'm trying to investigate whether there are errors within my BAM. After looking at the BAM header to see whether the read groups exists (using samtools, i.e.

samtools view -H file1.bam

Here, it appears that the header includes @RG tags. However, I suspect that the BAM contains reads that either lack the @RG tag or perhaps these read groups are not well-defined, i.e. it looks like some reads in the BAM may be missing an RG tag.

How do I explicitly check whether all reads have an RG tag, and whether these are well-defined?

According to the man page, running samtools stats --split RG <file1.bam> should produce summary statistics separated by read group. If it doesn't produce a list of ungrouped reads, counts/statistics can be compared to running without the --split RG argument.

Here's some example output from a BAM file with combined read groups:

$ samtools stats --split RG both.bam | grep '^SN'
SN  raw total sequences:    2392
SN  filtered sequences: 0
SN  sequences:  2392
SN  is sorted:  1
SN  1st fragments:  2392
SN  last fragments: 0
SN  reads mapped:   2341
SN  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
SN  reads unmapped: 51
SN  reads properly paired:  0   # proper-pair bit set
SN  reads paired:   0   # paired-end technology bit set
SN  reads duplicated:   0   # PCR or optical duplicate bit set
SN  reads MQ0:  0   # mapped and MQ=0
SN  reads QC failed:    0
SN  non-primary alignments: 0
SN  total length:   2292231 # ignores clipping
SN  bases mapped:   2255361 # ignores clipping
SN  bases mapped (cigar):   1538560 # more accurate
SN  bases trimmed:  0
SN  bases duplicated:   0
SN  mismatches: 254720  # from NM fields
SN  error rate: 1.655574e-01    # mismatches / bases mapped (cigar)
SN  average length: 958
SN  maximum length: 9078
SN  average quality:    20.0
SN  insert size average:    0.0
SN  insert size standard deviation: 0.0
SN  inward oriented pairs:  0
SN  outward oriented pairs: 0
SN  pairs with other orientation:   0
SN  pairs on different chromosomes: 0

The --split command has also produced additional .bamstat files, which can be looked at for the summary statistics for individual read groups. In this case I'm only interested in the "raw total sequences" line:

$ grep 'raw total sequences' *.bamstat
both.bam_minusOnly.bamstat:SN   raw total sequences:    1901
both.bam_plusOnly.bamstat:SN    raw total sequences:    491

There are 1901 sequences from the "minusOnly" group, and 491 sequences from the "plusOnly" group, and I know from the previous statistics that there are 2392 total sequences. 1901 + 491 = 2392, so I know there are no ungrouped reads in this file.


How to demultiplex pair-end fastq reads with barcode 2 in the identifier line? I have multiplexed pair-end fastq reads with dual barcodes. The issue is that one barcode is present in the header and one is present at the beginning of the read. I need a method to demultiplex this data, but in order to assign a read to an individual, both barcodes are required, as there is overlap between the barcodes. It seems there are packages available to demultiplex using header ID or in-line barcodes to demultiplex, but not both.

example reads:

@700819F:525:HT235BCXX:2:1101:1139:2144 1:N:0:ATCACGAT
CGAATTGCAGATTTTTTCTGAATAAAGCAGTGCAATAAAATTCCCCGCAAAAACACTTNANNNGNNNNNNNNNNNNNNNNNNNNANNNNNNGNTAATAAA
+
GGGGGIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#<###.####################.######<#.<<GGGI
@700819F:525:HT235BCXX:2:1101:1212:2172 1:N:0:ATCACGAT
AAGGATGCAGGGCATCTCCCTCAGGCTGCGCTCTATCGAAGTCATCCCAGAATTAGATTCCGACCACAGACCAGTCTTAGTCAAACTAGGACCCGAGTGT
+
GGGGGIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIGIIIIIIIIIIIIIIIIIIIIGIIIIIIGGIIIGIIIIIIIIIIIIIGGIGIGIIIGIGI
@700819F:525:HT235BCXX:2:1101:1110:2173 1:N:0:ATCACGAT
GGTTGTGCAGAAAGAGTTGCTGATAAACTTAGCCATGCAGAACAGAATTATGAGTTAGAAGTATGTATATATATACCAATCACTATATCAACCCATTACC
+
<G.G<GGIIIG.GA.GAGAAGG<AA<A<.<<<GA<.<<.G<.G<<A<GGAA....<G.G..<<.GA..<A.<GG<<<.<..<<.GG.A..G..<<<.<<G

One of the barcodes is the sequence found at the end of the header line (ATCACGAT). The second barcode is found at the beginning of the sequence (CGAAT in read 1). Both are required for assignment to an individual.

You could try something like that:

#!/bin/sh

fq_in="${1}"

bioawk -c fastx '{print substr($comment, 7, 15)}' ${fq_in} | uniq \
    > "barcode1.txt"
bioawk -c fastx '{print substr($seq, 1, 5)}' ${fq_in} | uniq \
    > "barcode2.txt"

for bc1 in $(cat "barcode1.txt")
do
    for bc2 in $(cat "barcode2.txt")
    do
        fq_out="${bc1}_${bc2}.fq"
        bioawk -c fastx -v bc1=${bc1} -v bc2=${bc2} \
            'substr($comment, 7, 15) == bc1 && substr($seq, 1, 5) == bc2 {print "@"$name" "$comment"\n"$seq"\n+\n"$qual}' \
            ${fq_in} \
            > ${fq_out}
    done
done

This might be slow. You can a little bit of time if you directly provide lists of barcodes for the for loops.

It could be more efficient to write a custom program that can write to multiple files. It happens that I wrote something that does this in python just yesterday, but only reading the inline barcode. Modifying it to also use the in-header barcode is likely simple.


Omics data: How to interpret heatmap and dendrogram output? How to interpret heat map and dendrogram output for biological data (omics) in words (when writing results and discussion)?

What should I consider (statistics behind?) and what is the best approach?

Here is one of my HM for proteomics data.

Script being used

col <- colorRampPalette(c("red","yellow","darkgreen"))(30)
png("HM1.png")
heatmap.2(as.matrix(df), Rowv = T, Colv = FALSE, dendrogram = "row", #scale = "row", col = col, density.info = "none", trace = "none", margins = c(7, 15) )
dev.off()

Heat Map output.

I would be greatful if you would give an example.

The dendrogram summarize the information of a group of values and sort them according to the similarity they have. It can be applied to both, samples and features.

The dendrogram allows to visualize features that are more similar together, usually revealing patterns that wouldn't have been seen otherwise.

In an article usually it is used something like: "we can see that X and Y are clustered together revealing that ... (they have something in common), while Z and A are clustered far"


CRISPR screening and systems biology I am not sure how to tackle this one, so I’ll explain the general idea I have in mind.

Given results of a CRISPR knockout screening experiment (like in this paper) with two types of assays: control vs. treated with 5 time points per each of the assays.

Is there a way to use the screening data to perform systems biology network analysis? Something which is related to signal transduction, transcriptional regulatory, or metabolic networks?

Maybe something which uses Bayesian inference? Or something similar to differentially expressed representation of the CRISPR screening? Alternatively, a GWAS analysis similar to this paper?

What do you call screening results?
Do you have a kind of array with protein or RNA data?
If so, you can use several open-source tools available on the internet, like Enrichr or network analyst to better understand the pathways affected by your treatment.


how do i change a variable with iterator under a for loop in python This question comes up in my bioinfo class.

I have a csv file containing gene names and the cluster they belong to. I want to run GESA on each individual cluster. So I would like to extract the list of gene names corresponding to each cluster. So I wrote code like this

import pandas as pd
gene_list = pd.read_csv('marker_genes.csv')
for i in range(12):
    "cluster_"+str(i) = gene_list.gene[gene_list.cluster==i]

I know this raises syntax error, so what is the correct way of doing it.
In the end, I would want cluster_0 to cluster_12 each being a list of gene names.

cluster_0 to cluster_12 are just strings. they cannot be the names of arrays or matrices. Strings can only be used in print statements and string processing statements.

The correct way to define an array or matrix is as follows: call the arrays as: cluster[0] to cluster[12]

then, assign:

  for i in range(12):
        cluster[i] = gene_list.gene[gene_list.cluster==i]

then whenever you need to print the names of these arrays, you use:  "cluster_"+str(i)


Question of Padj value (very high, infinite) got from dds After I got dds from my count matrix by DESeq(), I use results() method to get Padj and log2 fold changes for volcano plotting,
then I found some gene's Padj values are infinite and the P-value is very close to 0;
this makes those dots on my volcano plotting located very top to Y-aix and also confused me:

is this normal or not, such a p-value (very close to 0) and -logPadj is infinite;
how should I treat those data: are they very different expressed between samples and control or just abnormal data which should be taken off?

Here follows the table:

table 1 is dds results, -log10(padj) show inf;
table 2 is count matrix from featureCount
table 3 is size factor of samples
Please help me !

The fold changes for those are huge (on the order of 32x for the largest), so it's unsurprising for the adjusted p-values to be absurdly small.
There's no reason to treat them differently.


How do kmer counters determine which kmer is 'canonical'? When counting canonical kmers, ie kmers in which both the forward and reverse complement of a sequence are treated as identical, how do kmer counting programs decide which kmer to use as the canonical sequence? Do they all work the same way?

To investigate I made a string with GAGTGCGGAATACCACTCTT which contains all 16 possible 2mers. I then used kmc to figure out how they determine which kmer is used. Only the kmers in the filtered column below appeared. So, it looks like KMCs' 'canonical' kmers are the ones that first occur alphabetically.

╔════════════════╦═════╦════════════════════╦══════════╗
║ Possible Kmers ║ RCs ║ RC occurs earlier? ║ filtered ║
╠════════════════╬═════╬════════════════════╬══════════╣
║ TT             ║ AA  ║ YES                ║ TA       ║
║ TG             ║ CA  ║ YES                ║ GC       ║
║ TC             ║ GA  ║ YES                ║ GA       ║
║ TA             ║ TA  ║                    ║ CG       ║
║ GT             ║ AC  ║ YES                ║ CC       ║
║ GG             ║ CC  ║ YES                ║ CA       ║
║ GC             ║ GC  ║                    ║ AT       ║
║ GA             ║ TC  ║                    ║ AG       ║
║ CT             ║ AG  ║ YES                ║ AC       ║
║ CG             ║ CG  ║                    ║ AA       ║
║ CC             ║ GG  ║                    ║          ║
║ CA             ║ TG  ║                    ║          ║
║ AT             ║ AT  ║                    ║          ║
║ AG             ║ CT  ║                    ║          ║
║ AC             ║ GT  ║                    ║          ║
║ AA             ║ TT  ║                    ║          ║
╚════════════════╩═════╩════════════════════╩══════════╝

Do all kmer counting programs use the same canonical kmers, and if so do you have documentation explaining this? I wasn't able to find anything in the papers for jellyfish or kmc.

When a k-mer is identical to its reverse complement, both are canonical. Note that a canonical k-mer is a sequence, irrelevant to its position(s) in the input string. More precisely, give a string $s$, its canonical string is
$$
{\rm canonical}(s|h)=\left\{\begin{array}{ll}
s & \mbox{if $h(s)<h(\overline{s})$}\\
\overline{s} & \mbox{otherwise}\\
\end{array}\right.
$$
where $\overline{s}$ is the Watson-Crick reverse complement of $s$ and $h$ is an arbitrary string hash function. In practice, we most often use 2-bit encoding as $h$. Order under such a hash function is equivalent to the lexicographical order.


How to filter ribosomal RNA from scRNA-seq data I want to filter out ribosomal RNA from scRNA-seq data (downloaded from here). Is there a list of known ribosomal RNA?

The only solution I found is SortMeRNA, however it works with raw sequencing data afaik, while I already have a matrix with transcript counts for each cell. I searched for a comprehensive list of rRNAs but I didn't find any.

In the paper mentioned, we used the ScaleData function in Seurat to regress out the number of reads, Rn45s abundance, and percent ribosomal gene transcripts. Ribosomal genes were found with the regular expression ^Rp[sl][[:digit:]].

tiss <- ScaleData(object = tiss, vars.to.regress = c("nReads", "percent.ribo","Rn45s"))

Here's a fuller notebook, and we'll have a better organized repository soon.


RDAVIDWebService (R package) manual installation, or alternative tool Recently, I had to re-install all my R packages (R version 4.1.0 (2021-05-18) -- "Camp Pontanezen") and it turns out that "RDAVIDWebService" is not available anymore on Bioconductor:

package ‘RDAVIDWebService’ is not available for Bioconductor version '3.13'

Is there an another way to install it (couldn't find the files on the Bioconductor site to install the package manually) or is there another R package that allows to use the DAVID tool?

That package was deprecated in Bioconductor version 3.13 after several attempts to contact the package author/maintainer.
Your options include reverting back to Bioconductor version 3.12 or using one of the other options for GO/pathway enrichment analyses, e.g. topGO, enrichR, or clusterProfiler.


How did researchers derive the Ramachandran "validation" contours? I'm a beginner to structural biology and for fun, calculated the tortional angles of some 100, 000 proteins. Here is my Ramachandran plot:

When I went to look for "canonical" Ramachandran plots, I discovered most researchers overlay a "contour" plot of allowed torsional angles (see below). My two questions are:

How are the distributions (ie. countours) of protein secondary structures calculated/derived? (ie. How do I "cluster" these into beta-sheets, alpha helices?)
How can I derive them myself?

To expand on Matteo Ferla's comment as an answer, there are two approaches to getting these "allowed" contours:

From probabilities found in observed crystal structures. You could imagine plotting a Ramachandran plot for the whole PDB and colouring the points by assigned secondary structure from DSSP. The contours would then be regions that encompass e.g. 99% of alpha-helical residues, or 99% of all residues. In doing this you would notice for example that glycine can adopt more regions of the Ramachandran space that other amino acids. PROCHECK is software that scores a protein based partly on the likelihood of the phi/psi angles of each residue.
From reasoning about the chemistry of the secondary structure of proteins. Only certain phi/psi combinations allow the i -> i + 4 pattern in an alpha-helix or the alternating up/down sidechain pattern in a beta-strand.


Differential gene expression analysis of time series with replicates I have a dataset that has two groups, perturbed vs control. Each group has 3 replicates. Each replicate has 8 time points. How do we do Differential gene expression analysis to find significantly expressed genes? Any packages in R or web tools etc?

R packages like DESeq2, edgeR, and sleuth will help you perform differential gene expression analysis.
These are typically used for control-vs-treatment experimental designs.
So you could use them like this and perform a bunch of pair-wise comparisons between the time points (i.e. what are the differentially expressed genes between timepoint $t_0$ and timepoint $t_1$, $t_0$ and $t_2$, etc.).
But all of these methods make use of generalized linear models to perform the statistical inference on coefficients.
For the simple control-vs-treatment design, above, this looks something like
$$g \left(\mathbb{E} \left[Y_i | X \right] \right) = \mu_i + \beta_i X$$
where $g$ is the link function, $\mu_i$ is the mean expression of gene $i$, $\beta_i$ is the expression fold change estimate, and $X$ is the indicator variable saying if the sample is a treatment or control sample.
You can modify this equation by adding a time component to it.
$$g \left(\mathbb{E} \left[Y_i | X \right] \right) = \mu_i + \beta_i^{(group)} X + \beta_i^{(time)} t$$
If you specify the timepoint of the sample in the experimental design matrix for DESeq2, edgeR, or Sleuth, you can get information about the effect of the perturbation as well as the effect over time.
You can try more complicated designs, too, like a coefficient for the interaction between timepoint and perturbation status, but with 3 replicates in each, it's probably best to stick with simple, linear models.


Somatic tumor only variant calling? I'm evaluating possibilites for somatic tumor variant calling without paired-normal samples. I'm aware of the consequences without a normal sample.

All the popular variant calling tools such as Strelka, VarScan etc require a normal sample.

What are the some of the possibilites I have? I'm looking for good sensitivity in exchange for higher false positives.

There is a recent paper that attempts to do this:

  ISOWN: accurate somatic mutation identification in the absence of
  normal tissue controls.

  In this work, we describe the development, implementation, and
  validation of ISOWN, an accurate algorithm for predicting somatic
  mutations in cancer tissues in the absence of matching normal tissues.
  ISOWN is available as Open Source under Apache License 2.0 from
  https://github.com/ikalatskaya/ISOWN


Mapping RefSeq-ID or NCBI-ID (Entrez-ID) to EC number I annotated my bacterial genomes using the new NCBI Prokaryotic Genome Annotation Pipeline and now, I want to annotate EC-numbers.

In the master annotation file (.gff), I get the RefSeq-ID of every successful annotation, e.g. RefSeq:YP_805528.1.

When I visit the respective NCBI page, I see that it sometimes is mapped to an EC number.

I'd like to automate this process: Convert every RefSeq-ID to an EC number. How can I do this?

(I may already have discovered the first step: I can map the RefSeq-ID to the NCBI-ID (Entrez-ID) using this NCBI conversion database: YP_805528.1 -> 4419225. But I don't see how to get from these databases to EC.)

Are there other ways? As far as I know, most people blast their genes to already-mapped databases like SwissProt and ExPASy ENZYME. But I don't like this approach because technically, the NCBI annotation pipeline already did this. Why compute the same thing twice?

You can use Entrez Direct for this as follows:

efetch -db protein -id YP_805528.1 -format gpc -mode xml | xtract -insd Protein EC_number
YP_805528.1     6.3.2.13


Downloading multiple SRA files from several SRA accession IDs does not work I am trying to download multiple SRA files located in several SRA accessions. Some of my accession numbers are as follows:

SRR1002678
SRR1158898
SRR1158907
SRR1158950
SRR1181152
SRR1181155
SRR1181196
SRR1181204
SRR1181205
SRR1181208
SRR1269497
SRR2100227
SRR2100419
SRR2100830
SRR2101344
SRR671725
SRR671726
SRR671740
SRR671759
SRR671760
SRR671761
SRR671779
SRR671782
SRR671815
SRR671821
SRR671842
SRR671862
SRR671863
SRR671864
SRR671865
.............

I tried the following command:

for i in {SRR1002678..SRR1184123}; do wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR1*/$i; done

but this did not work.

I tried prefetch as well:

sort -u SRR.txt | parallel -j 4 "prefetch {}"

where SRR.txt is the list of my SRA IDs. However, with this command only a few files are downloaded as *.sra files, while with several other SRA IDs I get the following error:

2018-11-15T06:17:20 prefetch.2.8.2: 1) Downloading 'SRR1183962'...
2018-11-15T06:17:20 prefetch.2.8.2:  Downloading via https...
2018-11-15T06:20:21 prefetch.2.8.2: 1) 'SRR1183962' was downloaded successfully
2018-11-15T06:20:27 prefetch.2.8.2: 'SRR1183962' has 0 unresolved dependencies
2018-11-15T06:20:38 prefetch.2.8.2: 'SRR1183962' has no remote vdbcache

2018-11-15T06:18:59 prefetch.2.8.2: 1) Downloading 'SRR1184034'...
2018-11-15T06:18:59 prefetch.2.8.2:  Downloading via https...
2018-11-15T06:20:40 prefetch.2.8.2: 1) failed to download SRR1184034
2018-11-15T06:20:40 prefetch.2.8.2 sys: libs/kns/unix/syssock.c:606:KSocketTimedRead: timeout exhausted while reading file within network system module - mbedtls_ssl_read returned -76 ( NET - Reading information from the socket failed )
2018-11-15T06:20:40 prefetch.2.8.2 int: libs/kns/unix/syssock.c:606:KSocketTimedRead: timeout exhausted while reading file within network system module -

I do not have a range of SRA files to download but quite random SRA files to download. What am I doing wrong? Any help is highly appreciated.

First of all, the brace expansion you tried to use isn't valid. While bash can expand {1..3} to 1 2 3, it has no way of knowing how to expand {SRR1002678..SRR1184123}:

$ echo {1..3}
1 2 3
$ echo {SRR1002678..SRR1184123}
{SRR1002678..SRR1184123}

Perhaps you meant SRR{1002678..1184123}? But there's no need for this, you already have the accessions in a file and can read that. I see that the FTP server you are accessing has organized the data like this:

sra/sra-instant/reads/ByRun/sra/SRR/SRRNNN/SRRNNNXXXX/SRRNNNXXXX.sra

Where NNN are the 1st three digits of the accession. So, for example, the file for accession SRR1002678 is at:

ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR118/SRR1181152/SRR1181152.sra

With this in mind, you can write a little script that extracts the URL from each accession:

while read i; do
    n=$(echo "$i" | sed -E 's/(SRR...).*/\1/');
    wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/$n/$i/$i.sra
done < SRR.txt

The only trick is the sed command which simply looks for SRR and the three next characters and removes everything else on the line (it "captures" the found pattern and then replaces the line with what was captured (\1)).


Exon-exon junctions: compare experimental transcripts to reference annotation My aim is to parse an experimental transcript set (obtained by RNAseq) to check which splice junctions are already reported in a reference annotation and which ones are new.

I tried Regtools:

feed the experimental alignment file (BAM: aligned transcripts to reference genome) to Regtools junctions extract
use Regtools junctions annotate to compare with the reference annotation (GTF)

Exact commands I used for this are:

[..]/regtools/build/regtools junctions extract aln_sorted.bam -o regtools.bed - where the sorted BAM file has index in the same folder
[..]/regtools/build/regtools junctions annotate regtools.bed [..]/ref/genome.fa [..]/ref/annotation.gtf -o regtools_annotate_out

The alignments have been generated with Minimap2.

It did not behave as expected: this gives me as output a file with only junctions annotated as 'new' ($11==N), while other annotations should be possible according to the manual:

  DA - This exon-exon junction is present in the transcriptome provided by the user(GTF) The ends of this junction are hence known donor and known acceptor sites according to the GTF file.

  NDA - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The ends of this junction are
   known donor and known acceptor sites according to the GTF file.

  D - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The donor of this junction is a known donor but the acceptor is novel.

  A - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The acceptor of this junction is a known acceptor but the donor is novel.

  N - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The ends of this junction are hence not known donor/acceptor sites according to the GTF file.

My questions are:

am I doing something wrong?
do you have any experience with Regtools and can help debugging this behavior?
do you have other suggestion to achieve the output I am seeking?

You can use the script that comes with minimap2:

# install the k8 javascript shell
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` $HOME/bin/k8    # assuming $HOME/bin in your $PATH

# acquire the latest minimap2
git clone https://github.com/lh3/minimap2.git

# print novel/erroneous junctions
minimap2/misc/paftools.js junceval -e anno.gtf rna-seq.sam > output.txt

Some example output:

P       read1    4       chr2    227332588       227340301       [(227332588,227340291)]
P       read2    1       chr9    19376379        19376489        [(19376388,19376493)]
P       read2    3       chr9    19378506        19378707        [(19378514,19378707)]
N       read3    1       chr1    160647335       160647472
N       read3    2       chr1    160647581       160659672

Here column 4–6 gives junctions not in the annotation. Column 7, if not empty, gives the list of annotated junctions/introns that overlap the junction in read alignment. Note that a big fraction of outputted records are errors in read alignment.


What are the differences between "Site" and "CONECT" in a PDB file The CONECT section in a PDB file describes atomic connectivity.

The SITE section describes (among other characteristics), environments surrounding ligands present in the structure.

Now, how come not every atom in the CONECT section is described in the SITE section?

For example: 1a28.pdb

how come not every atom in the CONECT section is described in the SITE section?

Because it is not required and not needed.

CONECT is for listing bonds, so that molecular graphics programs know how to draw the structure.

SITE is for annotation of interesting sites. Each listed SITE has a corresponding REMARK 800 that "details its significance". A graphical program may provide a menu with a list of all interesting sites for easy navigation.

Note: PDB is a legacy file format. The primary file format used nowadays by the PDB (organization) is called PDBx/mmCIF. So when looking into such details it is often better to look directly into mmCIF files.


How to adapt the fgseaL function to perform rapidGSEA computation of gene ranks across 9 different phenotype labels? I wish to adapt the r language function fgseaL, https://github.com/ctlab/fgsea , to perform rapidGSEA, https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-016-1244-x , computation of inter-class  deviation per gene and the subsequent gene rank sorting operation on 9 different phenotype labels as illustrated in the diagram immediately below:

I thought of applying the R-language rank() function on the Expression Data Matrix D. If that is not correct, what sequence of R language commands should we
apply to the Expression Data Matrix D to calculate a key value sorted deviation measure  across 8 labeled human leukemia groups and a healthy labeled normal control group prior to running fgseaL?

I show below how fgseaL finds the correlation matrix between the R language variable , mat, which corresponds to the Expression Data Matrix D and the R language variable , labels , which is a vector of gene phenotype labels

    tmatSc <- scale(t(mat))
    labelsSc <- scale(labels)[, 1]

    minSize <- max(minSize, 1)

    pathwaysFiltered <- lapply(pathways, function(p) { as.vector(na.omit(fmatch(p, rownames(mat)))) })
    pathwaysSizes <- sapply(pathwaysFiltered, length)

    toKeep <- which(minSize <= pathwaysSizes & pathwaysSizes <= maxSize)
    m <- length(toKeep)

    if (m == 0) {
        return(data.table(pathway=character(),
                          pval=numeric(),
                          padj=numeric(),
                          ES=numeric(),
                          NES=numeric(),
                          nMoreExtreme=numeric(),
                          size=integer(),
                          leadingEdge=list()))
    }

    pathwaysFiltered <- pathwaysFiltered[toKeep]
    pathwaysSizes <- pathwaysSizes[toKeep]

    corRanks <- var(tmatSc, labelsSc)[,1]
    ranksOrder <- order(corRanks, decreasing=T)
    ranksOrderInv <- invPerm(ranksOrder)
    stats <- corRanks[ranksOrder]

    pathwaysReordered <- lapply(pathwaysFiltered, function(x) ranksOrderInv[x])

    gseaStatRes <- do.call(rbind,
                           lapply(pathwaysReordered, calcGseaStat,
                                  stats=stats,
                                  returnLeadingEdge=TRUE))

I found a problem with the algorithm shown immediately below.
correcttest <- data.frame(names = row.names(normal))
correcttest <- cbind(correcttest3, normal)
correcttest <- cbind(correcttest3, ALL3m)
rownames(correcttest) <- correcttest$names
correcttest$names <- NULL
correctlabelnormal <- rep(0:0, 73)
correctlabelALL3m <- rep(1:1, 122)
correctlabel <- as.vector(c(correctlabelnormal,correctlabelALL3m))
s <- apply(correcttest, 1, function(x) coef(lm(x~correctlabel))[2])
o <- rank(s)
o <- max(o) - o + 1
res <- fgseaL(df,o,correctlabel,nperm = 2000,minSize = 1, maxSize=50000)
empty data table (0 rows) of 8 columns:   pathway,pval,padj,ES,NES,nMoreExtreme,size

I found the binary phenotype labeled group fgseaL test results below looked satisfactory.
correcttest <- data.frame(names = row.names(normal))
correcttest <- cbind(correcttest3, normal)
correcttest <- cbind(correcttest3, ALL3m)
rownames(correcttest) <- correcttest$names
correcttest$names <- NULL
correctlabelnormal <- rep(0:0, 73)
correctlabelALL3m <- rep(1:1, 122)
correctlabel <- as.vector(c(correctlabelnormal,correctlabelALL3m))
fgseaL(df,correcttest,correctlabel,nperm = 2000,minSize = 1, maxSize=50000)
       pathway        pval        padj         ES       NES nMoreExtreme  size
1: Gene.Symbol 0.003940887 0.003940887 -0.2460126 -1.180009            3 45714
                                 leadingEdge
1: AKIRIN2,LRRC20,HSPA5,HSPA5,DTWD2,ZFYVE28,

Thank you for your consideration.

I finally found another answer to my question. Please read this great article in May 12 2017 BioMed Central (BMC) Bioinformatics article titled Ranking metrics in gene set enrichment analysis: do they matter?.

Also, please read this blog , Diving into Genetics and Genomics: Gene Set Enrichment Analysis (GSEA) explained.

After reading these two articles, my choice for the best rapidGSEA local ranking measure is Minimum Significant Difference (i.e., MSD), because it has the best overall false positive rate (i.e., FPR) for larger sample sizes.

Finally, it is important to realize that fgseaL's phenotype labeling can hypothetically be emulated by GSEA with either one of it's sixteen possible ranking metrics or a custom ranking metric.


How to extract metadata from NCBI's short read archive (SRA) for a few runs? I wish to extract metadata from a list of runs on NCBI's short read archive. For instance, I'd like to extract the library name ("HS0798") from the following run info:

https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=run_browser&run=SRR390728

I have to extract metadata for a few hundred run numbers. Is there a existing command-line tool (sratoolkit, python, perl, bash, etc) that will extract the metadata about the run?

I could grep it out of a webpage dump, but I notice there are XML schemas provided on the site, which to me indicates that there's some sort of XML API available somewhere.

Entrez Direct is perfect for this sort of thing. You can do something like this:

$ esearch -db sra -q 'SRR390728' \
  | efetch -format runinfo \
  | cut -f1,12 -d,
Run,LibraryName
SRR292241,HS0798
SRR390728,HS0798

The runinfo table format has 47 fields including experiment accession, project accessions, bioproject and biosample identifiers and such. You may want to explore that a little bit more if you are interested in downloading metadata in addition to the LibraryName.

Also, if you prefer to parse XML output, you can append -mode xml to the efetch command to retrieve the runinfo table in XML format. xtract tool that comes with the Entrez Direct package can then be used to convert that into a table.


How do I efficiently subset a very large line-based file? This has come up repeatedly recently: I have a very large text file (in the order of several GiB) and I need to perform line-based subsetting for around 10,000 lines. There exist solutions for specific scenarios (e.g. samtools view -s for randomly sampling BAM files) but sometimes my use-case doesn’t fit into these categories.

Unfortunately a naïve sed-based solution is extremely slow:

time sed -n -f <(awk -vOFS='' '{print $0, "p"}' line_numbers.txt) input_file > selected_lines.txt

Where line_numbers.txt is a file containing one line number per line.

Forget running this for 10,000 lines; it’s already grinding to a halt for a mere 1000.

How can I speed this up, ideally so that it scales only with the size of the input file, and has more or less constant runtime n the number of lines that I subset?

Turns out, simply keeping track of the next candidate line (after sorting the sample line numbers) fixes the performance issue, and most of the remaining slowness seems to be due to the overhead of actually reading the file so there’s not very much to improve.

Since I don’t know how how to do this in sed, and it’s not trivial in awk either, here’s a Perl script:

#!/usr/bin/env perl

use strict;
use warnings;

my $file = $ARGV[0];
my $lines_file = $ARGV[1];

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
chomp (my @lines = <$lines_fh>);
close $lines_fh;

@lines = sort {$a <=> $b} @lines;

open my $fh, '<', $file or die "Cannot read file $file";
my $line = 1;
my $next_line = 0;
while (<$fh>) {
    last if $next_line == scalar @lines;
    if ($line++ == $lines[$next_line]) {
        $next_line++;
        print;
    }
}
close $fh;

I’ve implemented a similar function in C++ for an R package, that's only slightly longer than the Perl script. It is ~3 times faster than the Perl script on my test file.


Program Quality control for data sequencing (nanopore) Hope all of you are fine. I am wondering if exists a program easy to use that can bring me Basic Statistics of Quality Control from raw data made by MinION device (nanopore).

I have read about FASTAQC and Porechop. My sequencing raw data was a result of a metagenomics analysis (16s and 18s). My main objective is only to filter the sequences and obtain a basic statistic about the lenght and the number of sequences ordered by the lenght. Thank you ladies and gentlemen

PD Sorry for my poor english. Day to day i'm trying to improve it.

A simple start would be to use MinIONQC, which is also available via bioconda under the name r-minionqc.


Preparing data for UpSet Plot in R As i have multiple comparison so I m trying this UpSet plot, but it needs data to be formatted in particular way, one of the ways is as a list, so now the issue is since between one comparison the no of genes coming out may not be same as the other one, so there is a difference in dimension so if i bind them as data frame then i get dimension error, so the best way to get through is make list so what i have tried so far

A <- read.csv("GMP_MONO_UP_DOWN_CDS_EPI.txt",header = TRUE)
B <- read.csv("CMP_GMP_UP_DOWN_EPI.txt",header = TRUE)
C <- read.csv("HSC_BLAST_UPDOWN_EPI.txt",header = TRUE,sep = '\t')

dim(A)
dim(B)
ONE <- A[c(1)]#############just taking out the gene column########
TWO <- B[c(1)]#############just taking out the gene column########
THREE <- B[c(1)]#############just taking out the gene column########
dim(ONE)
dim(TWO)
b1 <- apply(ONE,2,list)
b2 <- apply(TWO,2,list)

Now this is the code from their help page

upset(movies, nsets = 7, nintersects = 30, mb.ratio = c(0.5, 0.5),
      order.by = c("freq", "degree"), decreasing = c(TRUE,FALSE))

so the data "movies" class is data frame, there is not dimension issue as such as in my case if i want to intersect two comparison which may not have the same length, so how do i overcome the issue any suggestion or help would be highly appreciated. Update my question with one more data to add.

We can plot using lists as well, using fromList function, see example:

library(UpSetR)

myGeneSets <- list(
  set1 = c("gene1","gene2","gene3"),
  set2 = c("gene1","gene4","gene5", "gene7"),
  set3 = c("gene1","gene7")
  )

upset(fromList(myGeneSets))

  fromList

  A function to convert a list of named vectors to a data frame
  compatible with UpSetR.

fromList(myGeneSets)
#   set1 set2 set3
# 1    1    1    1
# 2    1    0    0
# 3    1    0    0
# 4    0    1    0
# 5    0    1    0
# 6    0    1    1


How do researchers account for heterozygosity in genome assembly? this is a very general question about diploid genome assembly. I am wondering how people deal with heterozygous region in a genome when assembling ?
Do they pick one of the haplotype to be in the referene and leave out the other ?
I am particularly thinking about situation similar to mammalian chromosome Y.
Do you put Y or X in the assembly ? An hybrid ?

Thanks

Generally a single allele is used. This may be arbitrary or match the major allele frequency of the species. For sex chromosomes, usually both/all are included, so for mammals both X and Y are part of the reference genome.


Gene expression Table to Expression Matrix converstion I have an RNAseq gene expression file (Count data) as follows, I need to conduct a Differential gene expression analysis between Patients, for that, I need to have this file as a Matrix of Rows as Genes and columns as counts (samples), how best can I get this as matrix file?
Thanks a lot for your help, Stay safe

Your question stated differently is as follows:

I have data in long format and need to change it to wide format. For that, I want to keep the gene_raw, use patient_id for the column names and gene_code for the row names. How can I do this easily?

Assuming you know how to read that into R as a data.frame that I'll call d:
library(tidyr)
library(dplyr)
reformatted = d %>% select(patient_id, gene_code, gene_raw) %>% spread(patient_id, gene_raw) %>% column_to_rownames(gene_code)

Or something along those lines.


Combat for multiplatform batch correction I have been looking ways to use data from various micro-array platform such as agilient, rosetta/merk. Some of the established method I came across is combat, but I'm not sure if it can be done for cross platform.

A similar question here where the answer is kind of can't be or shouldn't be done.

Any suggestion how to take account of batch effects across multiple platforms.

Usually with microarrays you want to make a case/control comparison, so I am going to assume that.

Data from different array platforms is generally difficult to compare:  each platform is measuring potentially a different part of the expression of the gene (different exons or different regions), each platform is likely to need different normalisation approaches, and have different dynamic range.

If you have case/control data from several different platforms, then the most robust thing to do is to analyse each platform separately, and produce a list of up and down-regulated genes according to each platform.

Then, it would be reasonably safe to combine the lists of genes in a 'meta-analysis' approach.  I would expect that in a lot of cases, 'up' and 'down' regulation would be correlated across platforms, but I would not necessarily expect very strong correlation in the magnitude of changes between platforms for individual genes.


What is the difference between fixed effects and random effects in the context of Linear-mixed models? The terms in LME i.e. fixed and random create confusion? What is the genesis that can distinguish between the two ?

I do not know much about statistics but I will try my best to explain.

First, random effects are defined as the factors (categories) in the population that we are not aware of (not observed), so we are randomly sampling levels of those factors when we sample the population.

Practically speaking, random effects can be found when there are hierarchical structures in the data. For example, we want to compare the relationship between students' SAT score and their college admission rates. You might just want to run a simple linear regression. But there might be some structure in the data that is hidden. For example, the students can be grouped into classes, and then schools. You can also group students by the social-economic status of their parents.

The idea is that students belonging to different categories may have different college admission rates even if they have the same SAT score. Therefore, it is important to consider some of these "random effects" when building your model.


mRNA cap identification by sequencing? If you sequence a positive-sense  RNA virus genome using NGS via Nanopore, Illumina or even the old Sanger method, is it possible to tell if the virus has 5' (five prime) cap at the 5' end of the RNA genome?

In short the answer is no, however the 'car-jumping method' described here it is possible, but it is about the enzymic MMLV reverse transcription in response to the capping structure. I assume that an alternative reverse transcriptase doesn't do 'car jumping' therefore any sequencing method will pick this up, but the difference between MMLV verse the heat-tolerant style reverse transcriptases. I would assume the secondary structure associated with viral capping prevents a low temperature reverse transcriptase, i.e. MMLV from correctly reproducing the vRNA into cDNA.
It is important to note that a cap, isn't just a virus thing, its mRNA eukaryotes per se. The cap is a 7-methylguanosine (m7G) linked to the 5′ end of the transcript by a 5′-5′ triphosphate bridge [m7G(5′)ppp(5′)N- or m7G(5′)ppp(5′)Nm in which Nm is a 2′-O methylated nucleoside].


Is there a way to measure cell line similarity using python? I have a list of cell lines that should be compared against another list of cell lines as pairs. I am trying to find a way in order for me to compare two cells such as MCF7    and 10964c cell lines. I am using python so, can anyone provide some steps so that I can do the comparison.

I also have their ontologies in an OBO format as I need them to do the comparison. if there is a library or a certain procedure that can solve this problem will be a great help.

Edit:

The comparison I'm trying to do is to measure the similarity in different ways, there are many measurement methods such as the distance in the ontological hierarchy between two cells by calculating the shortest path. You can also measure the number of common ancestors. But also, there are libraries that perform complicated (and more accurate) similarity measurements like JiangConrath, Jaccard, or Resnik which is popular. The problem with theses libraries is that they are not oriented to calculate cell line ontologies and may arequire annotation corpus for some similarity measures.

At its current form, your question is not clear, for example which aspects of these cell lines are your trying to compare?

Regarding your point on the ontologies, in Python you can use the pronto library to work with ontologies in OBO format.


Does Prokka or gff3 change the sequence / subject ID of fasta files? I used prokka to create a database for some fasta files but I noticed a strange difference between the prokka and fasta files. Normally the seqid for the gff3 output for prokka has seqid's that match the sequence id's of the fasta files. The files I just produced with prokka have a different seqid than my fasta files. My fasta files look like this

>NODE_108_length_645_cov_0.679537_ID_215

However the prokka output gff3 file for the same fasta file uses a seqid like this

gnl|X|JHIOGLGG_28

What would cause this difference? And do the last digits still correlate to the same seqid's even though the string before it has changed?

Prokka is an annotation tool. It takes your contigs file with fasta headers >NODE_108...and searchs for ORFs and their putative functions. Prokka produces submission ready output, the new fasta headers you see >JHIOGLGG_28...are locus tags and gene description from the annotation and are not related to contig node data. Please see the github information page for more info.


Using pysam with cython: htslib/kstring.h not found I'm trying to learn to use cython to speed up some code based on pysam. My issue is not strictly speaking about bioinformatics, but rather about building tools using a bioinformatics library. I hope this is still relevant in this site (I hesitated posting it as a github issue to pysam, alternatively, but I was not sure that this was a problem with pysam).

In my project's folder, the bam25prime/libbam25prime.pyx file has the following content:

# cython: language_level=3
from pysam.libchtslib cimport BAM_FREVERSE
from pysam.libcalignedsegment cimport AlignedSegment

cdef object ccollapse_ali(AlignedSegment ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    if (ali.flag & BAM_FREVERSE) != 0:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_end - 1),
            str(ali.reference_end),
            ".",
            ".",
            "-"])
    else:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_start),
            str(ali.reference_start + 1),
            ".",
            ".",
            "+"])

def collapse_ali(ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    return ccollapse_ali(ali)

In my setup.py, the setup function is called with this argument:
ext_modules = cythonize("bam25prime/libbam25prime.pyx").

When I try to do the build (running python3.6 setup.py build_ext), cython seems to miss some include paths, because I get the following error:

[1/1] Cythonizing bam25prime/libbam25prime.pyx
running build_ext
building 'bam25prime.libbam25prime' extension
gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -O3 -march=native -fomit-frame-pointer -fPIC -I/home/bli/.local/lib/python3.6/site-packages/pysam -I${HOME}/include -I/home/bli/include/python3.6m -I/home/bli/.local/include -c bam25prime/libbam25prime.c -o build/temp.linux-x86_64-3.6/bam25prime/libbam25prime.o
bam25prime/libbam25prime.c:497:28: fatal error: htslib/kstring.h: No such file or directory

If I manually re-run the gcc command, adding -I/home/bli/.local/lib/python3.6/site-packages/pysam/include/htslib, the header file gets found.

How do I get cython to correctly determine the path to the necessary header files?

I guess another way to formulate the question would be the following: How is cython supposed to determine what include flags to use?

Or maybe there is a problem with my pysam installation?

If I remember well, I installed it using pip3.6 --user install pysam. I have a ~/.pydistutils.cfg file containing the following settings

[install]
optimize=1
[build_ext]
include_dirs=${HOME}/include
library_dirs=${HOME}/lib
rpath=${HOME}/lib
user=1

Any advice appreciated.

Based on a suggestion by Devon Ryan, I searched how to set C_INCLUDE_PATH for cython, and found this issue. It refers to a (currently not working) feature present in cython's documentation, that should let libraries inform cython of the correct localisation for the header files in the module's setup:

ext_modules = cythonize(
    "bam25prime/libbam25prime.pyx",
    include_path=pysam_get_include())

As noted in the issue discussion, this seems to have no effect, but a workaround exists, consisting in using the Extension class to store the header localization information.

So the general answer is: Some python libraries (including pysam) provide a function that returns the paths to their include directories. This can be used to inform cython, but has to be done "manually".

My setup.py now looks as follows:

from setuptools import setup, find_packages
# If you have .pyx things to cythonize
from distutils.extension import Extension
from Cython.Build import cythonize
from pysam import get_include as pysam_get_include

name = "bam25prime"
__version__ = "0.1"

# https://github.com/cython/cython/blob/master/docs/src/reference/compilation.rst#configuring-the-c-build
extensions = [
    Extension(
        "bam25prime.libbam25prime", ["bam25prime/libbam25prime.pyx"],
        include_dirs=pysam_get_include()),
    ]

setup(
    name=name,
    version=__version__,
    # [description, author, email, licence, ...]
    packages=find_packages(),
    scripts=["bin/%s" % name],
    ext_modules = cythonize(extensions),
    # ext_modules = cythonize(
    #     "bam25prime/libbam25prime.pyx",
    #     include_path=pysam_get_include()),
    install_requires=["pysam", "pybedtools"])

Then, python3.6 setup.py build_ext succeeds.

And for whoever might be interested, here are the other files of the project:

bin/bam25prime is a link to bam25prime/bam25prime.py, which looks like this:

#!/usr/bin/env python3
# vim: set fileencoding=<utf-8> :
"""
When executed, this module collapses reads in a bam file to their 5-prime end
(reducing them to length 1) and writes the result in bed format.
"""

import argparse
from pybedtools import BedTool, Interval
from pysam import AlignmentFile
from .libbam25prime import collapse_ali

# pybedtools.Interval | pysam.AlignedSegment
# chrom               | reference_name
# start               | reference_start
# end                 | reference_end
# bool(int(strand))   | is_reverse
#
# Interval(reference_name, reference_start, reference_end, strand=strand)

def collapse_and_sort(alis):
    """
    Collapse elements from *ali* to their 5 prime end.

    The resulting :class:`pybedtool.BedTool` is sorted.

    :param alis: An iterable of :class:`pysam.AlignedSegment`
    :rtype: :class:`pybedtool.BedTool`

    :class:`pybedtool.BedTool` is an iterable:
    It can be iterated over in a for loop, or produce an iterator
    using :func:`iter` on it
    (see https://stackoverflow.com/a/28353158/1878788).
    """
    return BedTool(
        "\n".join(map(collapse_ali, alis)),
        from_string=True).sort(stream=True)

def main():
    """Main function of the program."""
    parser = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument(
        "-b", "--bamfile",
        required=True,
        help="Sorted and indexed bam file containing the aligned reads.")
    parser.add_argument(
        "-o", "--out_bed",
        required=True,
        help="Bed file in which the results should be written.")
    args = parser.parse_args()
    with AlignmentFile(args.bamfile, "rb") as bamfile:
        collapse_and_sort(bamfile.fetch()).saveas(args.out_bed)
    return 0

if __name__ == "__main__":
    exit(main())

bam25prime/libbam25prime.pyx:

#!/usr/bin/env python3
# vim: set fileencoding=<utf-8> :
# cython: language_level=3
"""
This library contains a cythonized version of a function to collapse aligned
segments to their 5' end.
"""

from pysam.libchtslib cimport BAM_FREVERSE
from pysam.libcalignedsegment cimport AlignedSegment

cdef object ccollapse_ali(AlignedSegment ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    if (ali.flag & BAM_FREVERSE) != 0:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_end - 1),
            str(ali.reference_end),
            ".",
            ".",
            "-"])
    else:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_start),
            str(ali.reference_start + 1),
            ".",
            ".",
            "+"])

def collapse_ali(ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    return ccollapse_ali(ali)

def main():
    """Example use case."""
    print("No examples to show here.")
    return 0

if __name__ == "__main__":
    import sys
    sys.exit(main())

(Any advice or suggestions welcome: I'd like to improve my python package-making skills.)

For those interested, to use cython with pybedtools I had to define the following Extension (where glob is glob.glob and /home/bli/src/bedtools2-2.27.1/ is the location of the source code of bedtools).

    Extension(
        "bam25prime.libcollapsebed", ["bam25prime/libcollapsebed.pyx"],
        include_dirs=glob("/home/bli/src/bedtools2-2.27.1/src/utils/*"),
        language="c++"),


How to downsample some of the samples in RNA-seq data? I have 40 samples and these are into two groups. I would like to perform a differential analysis between two groups. The library size of the samples is very low. But there are two samples (GroupA_12 and GroupB_2) in which the number of counts is too large compared to other samples. Please check the library sizes below:

I wanted to downsample these two cases and then would like to perform differential analysis. May I know how to downsample those two samples.

Any help is appreciated. thanq.

You have a few options:

Downsample the fastq files and rerun the entire analysis. You can do this with seqtk sample.
Downsample the BAM files, which you can do with samtools view -s.
Divide all of the counts in the counts files by some factor and round that to an integer.

I personally prefer option 2, since it's quick and doesn't usually have any side-effects unless you use two-pass alignment or did the quantification with something like salmon. Option 3 is always tempting, but I have to say I've never been convinced that this consistently produces similar enough results to options 1 or 2.


How can I subset WGS data to the level of WES variants? I would like to compare mutational signatures1 in patients from different studies, however some studies are based on exome seq (i.e. ~20,000 coding variants) and some are from whole genome seq (i.e. ~22,000 coding variants) - is there a way to 'downsample' WGS data to better reflect WES data and effectively 'ignore' the coordinates of those ~2000 coding variants in the VCF files?

1Alexandrov, L.B., Kim, J., Haradhvala, N.J., Huang, M.N.,
Ng, A.W.T., Wu, Y., Boot, A., Covington, K.R., Gordenin, D.A.,
Bergstrom, E.N. and Islam, S.A., 2020. The repertoire of mutational
signatures in human cancer. Nature, 578(7793), pp.94-101.

bcftools would be my choice, I'm sure bedtools could do the trick, too.
Something along this
bcftools view --regions-file

or  --targets, --regions-file might require a tabix index.


Convert sequence in MS doc to fasta or genbank How can I convert a sequence provided in a Microsoft doc file into a fasta or genbank format?

The CrazyDoc Python package can convert Genbank/Fasta/Snapgene/MSDoc sequences into Biopython records, which can be saved as Genbank or Fasta.
A web interface / demo is provided at EGF CUBA (Collection
of Useful Biological Apps): Convert Sequence Files.
Disclaimer: I'm the current maintainer of CrazyDoc.


Understanding the "relativeness" in RNA seq experiments Let's say we have two samples, sample 1 and sample 2, and we're running an RNA seq experiment. Consider genes A and B for both. For simplicity, let's say that when we run the RNA seq experiment, we get the following results:

gene A sample 1 count: 2

gene B sample 1 count: 4

gene A sample 2 count: 2000

gene B sample 2 count: 4000

Based on my understanding, RNA seq gives relative results. So I guess in this guess, how can you tell whether either one of the following are true:

1- You sequenced sample 2 a thousand fold more compared to sample 1

2- Sample 2 really has a thousand fold more expression in both genes A and B

We cannot be sure about this. For all we know, they could have been sequenced to the same depth and there is a biological difference that drives the 1000x (ridiculously large) fold change.
See above

The only thing you can say for sure is that gene B is expressed twice as much as gene A in samples 1 and 2. That's what "relative" means - relative expression of genes within a sample.


Error while calling bcftools mpileup - Failed to open -: unknown file type I have sequenced a bacterial genome with a GridIon from ONT. Basically what I want to check is whether or not trimming 50 bps at the beginning of the reads will improve alignment against the reference genome and ultimately the call of a consensus sequence.

For that what I have done so far is:

concatenated all the fastqs generated from the Gridion into a single (very large) fastq file: ref494_cat_all.fastq
aligned this concatenated fastq against the reference genome ref_genbank.fasta with minimap2: ~/Downloads/minimap2-2.17_x64-linux/minimap2 -ax map-ont ref_genbank/ref_genbank.fasta ref494_cat_all.fastq > ref494_aligned_trim50.sam
Performed the general samtools steps to covert sam to bam, index and sort

Then, in order to get the consensus sequence I first try to generate a .vcf.gz file, according to this tutorial

bcftools mpileup -f ref_genbank/ref_genbank.fasta ref494_aligned.sorted.bam | bcftools call -mv -Oz -o calls.vcf.gz

Immediatelly after calling the command above I got:

Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
[mpileup] 1 samples in 1 input files

Then, after a minute or two, I got:

Failed to open -: unknown file type

I'm pretty sure that there's no issue with .bam file, as I could visualise the alignment with IGV or tablet.

Does anyone know what the issue can be?

The - is usually used to mean standard input when reading data. So the message would suggest that the second bcftools command fails to read from stdin, so it fails to read the output of the first bcftools command:

$ echo "foo bar" | /ref/software/bcftools-1.8/bcftools call -mv -Oz -o calls.vcf.gz
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
Failed to open -: unknown file type

So there is something wrong with either your input file or with the output produced when running the first bcftools command on that input file. We can't really help more unless you share the file.


What are some of the major questions pursued by the field bioinformatics? I am an occasional user of tools in bioinformatics and I am fascinated by the possibility of uncovering deeper patterns in the biological realm using the help of computer algorithms. Yet I don't know much about the current state of this science. As I am considering a post-graduate degree in bioinformatics, I am interested to dig deeper into this subject and check whether the subject indeed aligns with my interest. So I start with the basic question stated in the title: What are some of the major questions pursued by the field bioinformatics?

Thanks in advance.

Welcome! Bioinformatics is just a term that describes an approach to doing biology. People who use ( DNA sequencing/genomics/protein informatics + programming ) to answer biological questions are after the same answers as the the people who use ( pipettes/microscopes/molecular biology + a lab notebook ).

My best suggestion is to learn about some biological concepts that are interesting to you, see how people have applied bioinformatics to that branch of biology, and think if there is something you might be interested in doing to work on that question. It is critical to have a strong understanding of the underlying biology to get to the core of new findings using bioinformatics.

A lot of the hardcore algorithm development going on today is in genome assembly with long + noisy DNA sequencing reads, aligning long+noisy reads to genomes and trying to figure out the differences between the maternal and paternal haplotypes, working on evolutionary questions with population-level genetic variation data, and figuring out how things are related to each other with phylogenetics, to name a few.

As far as a machine-learning type approach which you seem to be hinting at in your question, one very hot topic at the moment is single-cell sequencing. This is a set of technologies that allows us to get DNA, RNA, or protein sequencing from a population of single cells to look at heterogeneity in tissues, timepoints, et cetera.

Good luck! If you have some free time on your hand this website, ROSALIND, is a fun project Euler-type website where you can learn some biology while programming bioinformatics solutions.


Large-scale gRNA design for a CRISPR screen What are the best tools to design gRNAs in a high-throughput way for a CRISPR screen, e.g. targeting all protein-coding genes in a genome? I would like to take into account possible off-target effects, as well as to allow for flexibility in the PAM sequence, to have the possibility of using non-canonical PAMs.

I'm aware of CRISPRseek, but I'm afraid it might be quite slow to run for a genome-scale search.

Perhaps GenomeCRISPR, I've personally never used it, but it has a well-documented API and appears you could automate it to go through the whole genome.


What do the symbols mean in minimap2's gap cost equation? Heng Li's paper on minimap2 is understandable to me right up until the moment that it gives this definition of gap cost:

  2.2.1 Alignment with 2-piece affine gap cost

  Minimap2 performs DP-based global alignment between adjacent anchors in a chain. It uses a 2-piece affine gap cost (Gotoh, 1990):

  $$\gamma_a(l) = \min \{ q + |l| \cdot e, \tilde{q} + |l| \cdot \tilde{e} \}$$

Suddenly, out of nowhere, we've got five variables used without prior introduction: $l$, $q$, $e$, $\tilde{q}$, and $\tilde{e}$.

Heng doesn't state their meanings, and based on a quick skim and text search of the cited paper by Gotoh, that paper doesn't use these variable names anywhere either.

What does this equation mean?

Here is the translation:

$l$ = gap length
$q$ = initial gap penalty for shorter gaps
$e$ = short gap extension penalty per every additional gap
$\tilde{q}$ = initial gap penalty for longer gaps
$\tilde{e}$ = long gap extension penalty per every additional gap

To understand the first of the two expressions being minimized, see the Wikipedia section on affine gap penalties.

Notice the similar notation structure

$$A + B \cdot L$$

where

$A$ is the known gap penalty,
$B$ is the gap extension penalty, and
$L$ is the length of the gap

To understand the other variables with the tildes on top, see the text below Equation 2 (emphasis added to key explanation):

  On the condition that $e > \tilde{e}$,
  it applies cost
  $q + |l| \cdot e$
  to gaps shorter than
  $[(\tilde{q} - q)/(e - \tilde{e})]$
  and applies
  $\tilde{q} + |l| \cdot \tilde{e}$
  to longer gaps. This scheme helps to recover longer insertions and deletions (INDELs).

In other words, we make the assertion that we should penalize longer gaps differently than shorter gaps.

As a bonus, here's a visual understanding of why we use the minimization function on those two expressions using Python and some made-up (and probably not biologically useful) penalty terms. Note the requirement that $e > \tilde{e}$.

# Import libraries
import numpy as np
import matplotlib.pyplot as plt

# Setup
l = np.array(range(0, 15, 1))
y_short = eval('10 + 5*l') # Q=10, E=5
y_long = eval('20 + 3*l')  # Q=20, E=3

# Plot
plt.plot(l, y_short)
plt.plot(l, y_long)
plt.legend(["Short Gap Penalty", "Long Gap Penalty"])

Using our example and the expression $[(\tilde{q} - q)/(e - \tilde{e})]$, it will apply the first expression in the minimization function (the one without the tildes) when the gap length is shorter than 5.

$$
[(\tilde{q} - q)/(e - \tilde{e})]
=
[20 - 10)/(5 - 3)]
=
10 / 2
= 5
$$

We can see from our plot that the two lines intersect at 5, which we've verified through some arithmetic.


How can I extract the longest N isoforms per gene from a fasta file? The question of how one can extract the longest isoform per gene from a multi-fasta file has been answered expertly in this thread:
How can longest isoforms (per gene) be extracted from a FASTA file?

I wonder: how can the python script "filter.py" by Devon Ryan, that is marked as best answer, be modified so that it extracts the N longest isoform sequences per gene instead of only the single longest isoform?

Any help on this to a python novice is very much appreciated.

Here is an example of how the data looks (two genes with 5 isoformes each, and I need the 3 longest isoforms for each gene):

>Animal1|gene000001_1
TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK
>Animal1|gene000001_2
CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK
>Animal1|gene000001_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000001_4
HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF
>Animal1|gene000001_5
EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD
DPVVEIVLYLACSI
>Animal1|gene000002_1
TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK
>Animal1|gene000002_2
CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK
>Animal1|gene000002_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000002_4
LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF
ITSRSLR
>Animal1|gene000002_5
ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS

As output, I would like to have a fasta file like the following:

>Animal1|gene000001_5
EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD
DPVVEIVLYLACSI
>Animal1|gene000001_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000001_4
HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF
>Animal1|gene000002_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000002_4
LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF
ITSRSLR
>Animal1|gene000002_5
ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS

I know you are asking for a python solution. But maybe you are interested in other ways.

You can use samtools faidx to index your fasta file. The resulting file contains the information about the length of every sequence. We can sort by length and extract the first n transcripts per gene.

1. Index fasta

samtools faidx input.fa

2. Extract IDs of the n longest transcripts per gene

sed 's/_.\+//g' input.fa.fai | sort -u |
    parallel -k 'grep {} -m 3 <(sort -k2,2 -r --sort g input.fa.fai) |
        cut -f1' > longestID.txt

The sed command replaces everything from _ onwards in each line, leaving behind only the gene name. sort -u is used to get a list of distinct gene names. For each of theses gene names we grep for the first 3 matches in the descending sorted by length index file. As we just want the sequence ids, we select only the first column with cut. Adopt the -m parameter if you want more or fewer transcripts.

3. Get sequences

samtools faidx input.fa -r longestID.txt > output.fa


Does GISTIC (v 2.0) estimate amplified/deleted probabilities on a single sample basis? Does GISTIC 2.0 estimate the background model:

G = -log(Probability | Background)

by permuting within the sample or across all samples in the set?

The paper describes the probabilistic scoring method based on permutations, but I could not understand if this permutation is performed within the sample only. The documentation page seems to suggest across samples, but this would mean that different sized sets might lead to different outcomes on the sample sample.

Basically, does it matter if the set is constituted by, for example, 10 samples of the same tissue of origin or ~1000 from multiple tissues?

Thanks,

There doesn't seem to be much differences between GISTIC 1.0 and 2.0 as it says:

  As with GISTIC 1.0, we obtain P-values for each marker by comparing the score at each locus to a background score distribution generated by random permutation of the marker locations in each sample

But on the supplementary material of the GISTIC 1.0 there is a more detailed explanation of the method. See the "Stage 2" section:

  Second, we compare these G-scores to the distribution of scores expected if only random aberrations were observed. This distribution can be determined by rescoring the genome after permuting marker locations within each sample; we instead derive a semiexact estimate.

Furthermore in another section ("Stage 2: Aggregation of Data from Different Tumors to Differentiate Between Driver and Passenger Aberrations") it says:

  To determine which of the aberrations identified in Stage 1 are likely to represent driver events, we aggregate the data from all tumors used in the analysis to generate summary scores for amplifications, deletions, and LOH. The statistical significance of each score is determined by comparison to the distribution of scores obtained by all permutations of the data (using a semiexact approximation), with correction for multiple hypothesis testing.

However, the relevant section (of the supplementary materials) seems to be "Null Hypothesis Generation: An Analytic Derivation of the Null Distribution", where it describes the semiexact approximation used.

In a supplementary file it describes it as "Generate all permutations of SNP labels within each sample to simulate datasets with random aberrations"

Conclusion:

It seems that the background probabilities are calculated within sample.

I can't say if it matters how the sets are constituted but I would say that the broader the sets are, the better estimation of the structural variations it will perform but the same estimation of the background will be done.

Ultimately you can check the code of the program, or test with 10 samples and replace one of them to see if the results are changing accordingly.


How do PCR duplicates arise and why is it important to remove them for NGS analysis? I am trying to understand PCR duplicates in NGS analyses (actually whole-genome). I searched, and the best answer I found is in this blog.

However I don't understand if I understood how PCR duplicates arise  correctly because I cannot see the problem of having them in the downstream analysis - aside from computational problems, i.e. unnecessary redundancy.

If I understood correctly, PCR duplicates arise during library preparation when PCR amplifies the fragments with adapters. In this case, if you have duplicate fragments, you'll end up amplifying some of them twice or more.

However, when you are doing the mapping, they should fall on the same region, most likely decreasing the quality of the mapping (since they increase the consensus on a specific sequence, which could have been subject to sequencing errors). But aside from that, you should have the same mapping with respect to the mapping where duplicates are removed, although with decreased quality.

Is the quality problem the actual reason for removing PCR duplicates or is there anything that I'm missing?

In any scenario where depth of coverage is an important factor, PCR duplicates erroneously inflate the coverage and, if not removed, can give the illusion of high confidence when it is not really there.

For example, consider the following hypothetical scenario.

                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
               ..........................C.
                ............................
                      .............T..............
                       ..................C.........
                           ............................
                               ............................
                                   ............................
                                    .....C......................
                                       ..C.........................

The marked position has a coverage of 9 (9 overlapping reads), 4 of which suggest the presence of an alternate allele. Since these 4 reads map to different positions, they are independent observations supporting the same putative single nucleotide variant (SNV). The T in the third read can safely be ignored as a sequencing error, since it only appears once.

Now consider the following hypothetical scenario.

                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
                ............................
                      .............T..............
                       ..................C.........
                       ..................C.........
                       ..................C.........
                       ..................C.........
                               ............................
                               ............................
                               ............................

In this case, the marked position also has a coverage of 9, but 7 of those reads appear to be PCR duplicates. If we remove the duplicates, the marked position has only a coverage of 4, and only one C at the marked position, which is not enough confidence to call the SNV.


R package development: How does one automatically install Bioconductor packages upon package installation? I have an R package on github which uses multiple Bioconductor dependencies, 'myPackage'

If I include CRAN packages in the DESCRIPTION via Depends:, the packages will automatically install upon installation via devtools, i.e.

devtools::install_github('repoName/myPackage')

This is discussed in Section 1.1.3 Package Dependencies, in Writing R Extensions

Is there a way to streamline this such that packages from Bioconductor are automatically installed as well?

Normally, users install Bioconductor packages via BiocLite, e.g.

source("http://www.bioconductor.org/biocLite.R")
biocLite("edgeR")

As suggested, here’s an example showing the relevant lines from a DESCRIPTION file from a CRAN/GitHub hosted project that has Bioconductor dependencies (truncated):

Depends:
    R (>= 3.3.0)
biocViews:
Imports:
    methods,
    snpStats,
    dplyr

The relevant bit is the empty biocViews: declaration, which allows the Bioconductor dependency {snpStats} to be automatically installed.


How to download dbSNP database? I am bench-marking the performance of ADAM/Avocado on NGS whole-genome sequencing datasets. I am trying to repeat the tests reported in Frank Nofthat's (author of ADAM and Avocado) thesis, however, I am running into an issue acquiring the dbSNP database to use for Base Quality Score Recalibration.

The Base Quality Score Recalibration command (page 128-129):

adam-submit -- transformAlignments aln.adam bqsr.adam -recalibrate_base_qualities -known_snps dbsnp.adam

requires dbsnp.adam, which is described as "the dnSNP database" (converted to ADAM format) with reference to the NCBI database of genetic variation. This "database" can be found at this FTP site however I do not know which parts of these files to download and convert to ADAM in order to perform this step of the transformation pipeline.

Also, in order to use this genetic variation databaes do I need to make sure that the reference genome of "known variants" is the same genome that the sequencing data was aligned to?

Thank you.

In the dbSNP FTP site you linked, you need to go into organisms and select your organism of interest (human obviously). You can then select which release of the SNP database you want (eg. 150 vs 151) [the modified dates hint at when they were released] and which genome build you want (eg. hg37 vs hg38).

Below I am looking at the .vcf files. I assume this what you want because they can be converted to ADAM format via the adam vcf2adam command.

For example human (newest genome and newest dbSNP release):
ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/

Your next decision is whether or not you want to use ALL the annotated SNPs in the database or just the common ones. For the former you choose the All_20180418 file and for the later the common_all_20180418 file. I'm pretty sure the convention here is that the date in the filename is when it was uploaded/updated.

In your case I'm not sure what makes the most sense. Firstly, are you trying to replicate the initial conditions of when that software was published or do you want simply to benchmark the method with the best available data. Besides that, I imagine using all them is computationally much more taxing and if you are just coming up with general trend then using the common ones seems viable.