Toptable error, wont recognize condition I am getting a rather strange error from topTable.    When I run my code I get Error in fit$coefficients[, coef] : subscript out of bounds from topTable as if it is not recognizing the condition I am asking it to return the top ten DEGs of. I have provided a bit more context and some descriptions of my variables below.  Let me know if you need any more information or if anything is unclear.
Variables:
dge is a DGEList object of RNA seq data.
dge$counts is the raw counts.
dge$samples$dge_meta is the metadata and log2 transformed counts.
> colnames(dge$samples$dge_meta)
 [1] "Lab.sample.ID"
 [2] "Accession.Number"
 [3] "Alias.used.in.some.figures"
 [4] "Developmental.stage"
 [5] "Infected"
 [6] "Batch"
 [7] "Trimmed"
 [8] "Number.of.reads.that.pass.Illumina.filter"
 [9] "Total.number.of.reads.mapped"
[10] "%.of.total.reads.mapped"
[11] "Reads.mapped.to.T..cruzi.Esmeraldo.haplotype"
[12] "Reads.mapped.to.hg19"
[13] "%.of.mapped.reads.belonging.to.T..cruzi.Esmeraldo.haplotype"
[14] "%.of.mapped.reads.belonging.to.hg19"
[15] "sample"
[16] "hpi"

Lastly keep_it is a binary vector used to filter our zero count genes.
Code:
mm <- model.matrix(~~hpi*Infected, data = dge$samples$dge_meta) #making our design matrix

more_than_1cpm_model <- voom(dge$counts[keep_it,], mm, plot = TRUE) #this one is filtered by the 25% greater than 1 cpm

cutoff <-  1e-05

empBayesfit <- eBayes(fit_gt1cpm)

topTable(empBayesfit,
         coef = c("Infected"),
         p.value = cutoff, number = 10)

Error message:
Error in fit$coefficients[, coef] : subscript out of bounds

To understand the context better I ran rlang::last_error and rlang::last_trace as well.
> rlang::last_error()
<error/dplyr_error>
arrange() failed at implicit mutate() step.
* Problem with `mutate()` input `..1`.
x Input `..1` must be a vector, not a `standardGeneric` object.
i Input `..1` is `samples`.
Backtrace:
Run `rlang::last_trace()` to see the full context.
> rlang::last_trace()
<error/dplyr_error>
arrange() failed at implicit mutate() step.
* Problem with `mutate()` input `..1`.
x Input `..1` must be a vector, not a `standardGeneric` object.
i Input `..1` is `samples`.
Backtrace:
     x
  1. +-dplyr::arrange(x, samples)
  2. +-dplyr:::arrange.data.frame(x, samples)
  3. | \-dplyr:::arrange_rows(.data, dots)
  4. |   +-base::withCallingHandlers(...)
  5. |   +-dplyr::transmute(new_data_frame(.data), !!!quosures)
  6. |   \-dplyr:::transmute.data.frame(new_data_frame(.data), !!!quosures)
  7. |     +-dplyr::mutate(.data, ..., .keep = "none")
  8. |     \-dplyr:::mutate.data.frame(.data, ..., .keep = "none")
  9. |       \-dplyr:::mutate_cols(.data, ...)
 10. |         +-base::withCallingHandlers(...)
 11. |         \-mask$eval_all_mutate(quo)
 12. +-dplyr:::abort_glue(...)
 13. | +-rlang::exec(abort, class = class, !!!data)
 14. | \-(function (message = NULL, class = NULL, ..., trace = NULL, parent = NULL, ...
 15. |   \-rlang:::signal_abort(cnd)
 16. |     \-base::signalCondition(cnd)
 17. +-(function (e) ...
 18. | \-rlang::abort(...)
 19. |   \-rlang:::signal_abort(cnd)
 20. |     \-base::signalCondition(cnd)
 21. \-(function (cnd) ...

The error is telling you that there isn't any Infected column on the model.matrix.
Check the column names of mm. Also not sure that having two ~ will work well.


Solutions to reference bias issue I've recently learned about the reference bias issue - inability to properly map NGS reads in certain genomic regions caused by the fact that our standard genomic references are linear and do not fully represent the genomic variation present in the population. According to what I read this has multiple consequences: missing structural variants or variants in highly variable regions like HLA for humans.
Now I wonder what is the best way to address this. Of course, "best" is context-dependent but I don't have any specific context so general pros and cons will be highly appreciated. I know about two approaches but probably there are more. The first one is an alt-aware mapping that's used in DRAGEN by Illumina. As far as I understand, they just keep multiple linear alternatives for some regions and try doing a "smarter" mapping for those regions. The second approach that I've found is using graph references where the reference is not linear anymore but rather a graph and mapping is taking into account this graph structure.
Thanks in advance.

Answer from @throckmorton, converted from comment:
Here is a recent approach to this problem, using multiple linear references:
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02229-3


calling diploid SNVs from long reads I'd like to call diploid SNV variants from long-reads data (~80SMRTcells PacBio).

I have generated a draft reference genome for an indivudual from a heterozygous (~4%) species (Canu+Haplomerger2).

I can use this reference for some applications:

RNAseq mapping;
functional annotation;
duplicate detection;
synteny analysis

but not others:

analysis of compound heterozygotes;
allele-specific expression;
variant linkage

This is mainly why I need phasing. FYI, I tried Falcon + FalconUnzip too, but the daligner phase was just taking too much time/resources and I need to find a strategy to make it work in reasonable time (working on it).

Meanwhile, for phasing, It was suggested to use HapCUT2 - from Canu's FAQ

  (1) "Avoid collapsing the genome" in Canu to get raw diploid assembly as complete & separate as possible (2) HaploMerger2 to get high-quality reference and alternative haploid assemblies (3) HapCUT2 or other phasing tools to get the high-quality haplotype assembly based on the reference haploid assembly.

However, HapCUT2 requires the input of a diploid SNV file. Which could be done if I had short-reads data, but yet haven't found a method for doing so with in my hands just long-reads data.

Strategies I tried to get the SNVs:

marginPhase : method paper

but I had a issue in running it

Clairvoyante : method paper

but I currently don't have access to suitable GPU machines to build the model - and I can't use the models provided by the authors because my organism is too different from human

Any input would be greatly appreciated!

You could try aligning your reads to the draft reference genome with for example minimap2 and calling variants with freebayes. It appears there is a protocol for long reads: https://github.com/ekg/freebayes#re-genotyping-known-variants-and-calling-long-haplotypes

Disclaimer: I have not run freebayes myself, and it may be that this is a terrible idea. I would expect that calling SNVs (Single-nucleotide variants) from PacBio reads is very challenging.  However, with enough coverage, I think it could be possible.


What software/tool can I use to visualise a multiple sequence alignment in the xmfa format as produced by command-line tool Mauve? What software/tool can I use to visualise a multiple sequence alignment in the xmfa format as produced by command-line tool Mauve?

I am particularly interested in seeing the single nucleotide polymorphisms on the screen.

I am doing a multiple sequence alignment of about a dozen whole genomes, all in a single contig.

Link to Mauve

Typically I use the Mauve viewer to visualize the alignments.  Are there things you find lacking in that viewer?

There are multiple threads about converting xmfa to other formats


Can Maven (mass spectrometry software) read mzML files? I did not find a documentation. Can Maven (mass spectrometry software) read mzML files?

No, Maven still can read only mzXML files.


How can I extract normalized read count values from DESeq2 results? The results obtained by running the results command from DESeq2 contain a "baseMean" column, which I assume is the mean across samples of the normalized counts for a given gene.

How can I access the normalized counts proper?

I tried the following (continuing with the example used here):

> dds <- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno_treat)
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))

Here is what I have for the first gene:

> res["WBGene00000001",]$baseMean
[1] 181.7862
> mean(assays(dds)$mu["WBGene00000001",])
[1] 231.4634
> mean(assays(dds)$counts["WBGene00000001",])
[1] 232.0833

assays(dds)$counts corresponds to the raw counts. assays(dds)$mu seems to be a transformation of these counts approximately preserving their mean, but this mean is very different from the "baseMean" value, so these are likely not the normalized values.

The normalized counts themselves can be accessed with counts(dds, normalized=T).

Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. Given your previous question, we can see that geno_treat has a bunch of levels, which means that expanded models are not in use. In such cases, the baseMean should be the mean of the base factor in geno_treat.


How to generate the json files for Cromwell workflow execution? I am starting to use Cromwell WDL scripts for running GATK workflows on a slurm cluster. When using the paired-fastq-to-unmapped-bam script I need to define the fastq sepcifics in the appropriate json file. This needs to be written separately for each pair of fastqs. My question is: Is there a way of generating the json files on a batch of fastq-files by scripting? Of course, I could start writing scripts myself, but since I assume that GATK is used in a batch environment, I have the hope that something already exists.

Unfortunately generating inputs for Cromwell is not something that is handled by Cromwell itself or by GATK, and generally needs to be scripted to fit your custom needs.

Despite this being a bit arduous, there is a very good reason for this. WDL itself is a highly polymorphic workflow language, where inputs can be vastly different between each workflow. Cromwell ONLY understands how to parse the workflow and the inputs, but does not know the context of what it is running (ie this is a paired mapping job), nor does GATK know that it is running within a WDL.

Additionally, I've encountered many different cases where there are not just 2 fastq files but, due to older sequencing technologies, 10 or more for the same sample. Unless GATK/Cromwell knows how your specific samples are arranged there is no easy way for them to generate your inputs.

There are tools that would possibly help you simplify scripting better. WOMTool developed by the Broad has an inputs command which will take an  arbitrary workflow and generate a JSON template of the inputs, which you can then use in your own scripts.

As an example, imagine you have this WDL workflow in a file called my.wdl:

version 1.0
workflow foo {
  input {
    Array[Pair[File,File]] fastqFiles
  }
}

When you run WOMtool's inputs command, it produces a template JSON listing all inputs and specifying the variable types expected for each:

$ java -jar womtool.jar inputs my.wdl
{
  "foo.fastqFiles": "Array[Pair[File, File]]"
}

You can use this template to produce actual JSONs for each of your samples.

Note that you can also get the template written to a file as follows:

$ java -jar womtool.jar inputs my.wdl > my.inputs.json


Translating a genome sequence into possible frames I have a sequence below from MYC gene about which I need to translate the sequence in all possible frames AND identify (for each frame) which codons are actually used in MYC

AAAGCAATAATACAATTTAAAACCTGGGTCTCTAGAGGTGTTAGGACGTG
GTGTTGGGTAGGCGCAGGCAGGGGAAAAGGGAGGCGAGGATGTGTCCGAT
TCTCCTGGAATCGTTGACTTGGAAAAACCAGGGCGAATCTCCGCACCCAG
CCCTGACTCCCCTGCCGCGGCCGCCCTCGGGTGTCCTCGCGCCCGAGATG
CGGAGGAACTGCGAGGAGCGGGGCTCTGGGCGGTTCCAGAACAGCTGCTA
CCCTTGGTGGGGTGGCTCCGGGGGAGGTATCGCAGCGGGGTCTCTGGCGC
AGTTGCATCTCCGTATTGAGTGCGAAGGGAGGTGCCCCTATTATTATTTG
ACACCCCCCTTGTATTTATGGAGGGGTGTTAAAGCCCGCGGCTGAGCTCG
CCACTCCAGCCGGCGAGAGAAAGAAGAAAAGCTGGCAAAAGGAGTGTTGG
ACGGGGGCGGTACTGGGGGTGGGGACGGGGGCGGTGGAGAGGGAAGGTTG
GGAGGGGCTGCGGTGCCGGCGGGGGTAGGAGAGCGGCTAGGGCGCGAGTG
GGAACAGCCGCAGCGGAGGGGCCCCGGCGCGGAGCGGGGTTCACGCAGCC
GCTAGCGCCCAGGCGCCTCTCGCCTTCTCCTTCAGGTGGCGCAAAACTTT
GTGCCTTGGATTTTGGCAAATTGTTTTCCTCACCGCCACCTCCCGCGGCT
TCTTAAGGGCGCCAGGGCCGATTTCGATTCCTCTGCCGCTGCGGGGCCGA
CTCCCGGGCTTTGCGCTCCGGGCTCCCGGGGGAGCGGGGGCTCGGCGGGC
ACCAAGCCGCTGGTTCACTAAGTGCGTCTCCGAGATAGCAGGGGACTGTC
CAAAGGGGGTGAAAGGGTGCTCCCTTTATTCCCCCACCAAGACCACCCAG
CCGCTTTAGGGGATAGCTCTGCAAGGGGAGAGGTTCGGGACTGTGGCGCG
CACTGCGCGCTGCGCCAGGTTTCCGCACCAAGACCCCTTTAACTCAAGAC
TGCCTCCCGCTTTGTGTGCCCCGCTCCAGCAGCCTCCCGCGACGATGCCC
CTCAACGTTAGCTTCACCAACAGGAACTATGACCTCGACTACGACTCGGT
GCAGCCGTATTTCTACTGCGACGAGGAGGAGAACTTCTACCAGCAGCAGC
AGCAGAGCGAGCTGCAGCCCCCGGCGCCCAGCGAGGATATCTGGAAGAAA
TTCGAGCTGCTGCCCACCCCGCCCCTGTCCCCTAGCCGCCGCTCCGGGCT
CTGCTCGCCCTCCTACGTTGCGGTCACACCCTTCTCCCTTCGGGGAGACA

Where should I look for all 4 pieces of information easily? Please guide.

The 6-frame translation  results are below you will need to do the rest. It is translated 5' to 3' but appears to contain two splice sites. You will need to use blast to resolve the start/stop of the splicing.

All 6-codon translation is available the Expasy server here. You will need to blast to resolve the splicing.

5'3' Frame 1
KAIIQFKTWVSRGVRTWCWVGAGRGKGRRGCVRFSWNR-LGKTRANLRTQP-LPCRGRPRVSSRPRCGGTARSGALGGSRTAATLGGVAPGEVSQRGLWRSCISVLSAKGGAPIII-HPPCIYGGVLKPAAELATPAGERKKKSWQKECWTGAVLGVGTGAVEREGWEGLRCRRG-ESG-GASGNSRSGGAPARSGVHAAASAQAPLAFSFRWRKTLCLGFWQIVFLTATSRGFLRAPGPISIPLPLRGRLPGFALRAPGGAGARRAPSRWFTKCVSEIAGDCPKGVKGCSLYSPTKTTQPL-GIALQGERFGTVARTARCARFPHQDPFNSRLPPALCAPLQQPPATMPLNVSFTNRNYDLDYDSVQPYFYCDEEENFYQQQQQSELQPPAPSEDIWKKFELLPTPPLSPSRRSGLCSPSYVAVTPFSLRGD
5'3' Frame 2
KQ-YNLKPGSLEVLGRGVG-AQAGEKGGEDVSDSPGIVDLEKPGRISAPSPDSPAAAALGCPRARDAEELRGAGLWAVPEQLLPLVGWLRGRYRSGVSGAVASPY-VRREVPLLLFDTPLVFMEGC-SPRLSSPLQPARERRKAGKRSVGRGRYWGWGRGRWRGKVGRGCGAGGGRRAARARVGTAAAEGPRRGAGFTQPLAPRRLSPSPSGGAKLCALDFGKLFSSPPPPAAS-GRQGRFRFLCRCGADSRALRSGLPGERGLGGHQAAGSLSASPR-QGTVQRG-KGAPFIPPPRPPSRFRG-LCKGRGSGLWRALRAAPGFRTKTPLTQDCLPLCVPRSSSLPRRCPSTLASPTGTMTSTTTRCSRISTATRRRTSTSSSSRASCSPRRPARISGRNSSCCPPRPCPLAAAPGSARPPTLRSHPSPFGET
5'3' Frame 3
SNNTI-NLGL-RC-DVVLGRRRQGKREARMCPILLESLTWKNQGESPHPALTPLPRPPSGVLAPEMRRNCEERGSGRFQNSCYPWWGGSGGGIAAGSLAQLHLRIECEGRCPYYYLTPPLYLWRGVKARG-ARHSSRREKEEKLAKGVLDGGGTGGGDGGGGEGRLGGAAVPAGVGERLGREWEQPQRRGPGAERGSRSR-RPGASRLLLQVAQNFVPWILANCFPHRHLPRLLKGARADFDSSAAAGPTPGLCAPGSRGSGGSAGTKPLVH-VRLRDSRGLSKGGERVLPLFPHQDHPAALGDSSARGEVRDCGAHCALRQVSAPRPL-LKTASRFVCPAPAASRDDAPQR-LHQQEL-PRLRLGAAVFLLRRGGELLPAAAAERAAAPGAQRGYLEEIRAAAHPAPVP-PPLRALLALLRCGHTLLPSGR
3'5' Frame 1
CLPEGRRV-PQRRRASRARSGG-GTGAGWAAARISSRYPRWAPGAAARSAAAAGRSSPPRRSRNTAAPSRSRGHSSCW-S-R-GASSREAAGAGHTKREAVLS-RGLGAETWRSAQCAPQSRTSPLAELSPKAAGWSWWGNKGSTLSPPLDSPLLSRRRT--TSGLVPAEPPLPREPGAQSPGVGPAAAEESKSALAPLRSRGRWR-GKQFAKIQGTKFCAT-RRRREAPGR-RLREPRSAPGPLRCGCSHSRPSRSPTPAGTAAPPNLPSPPPPSPPPVPPPSNTPFASFSSFSRRLEWRAQPRALTPLHKYKGGVK---GHLPSHSIRRCNCARDPAAIPPPEPPHQG-QLFWNRPEPRSSQFLRISGARTPEGGRGRGVRAGCGDSPWFFQVNDSRRIGHILASLFPCLRLPNTTS-HL-RPRF-IVLLL
3'5' Frame 2
VSPKGEGCDRNVGGRAEPGAAARGQGRGGQQLEFLPDILAGRRGLQLALLLLLVEVLLLVAVEIRLHRVVVEVIVPVGEANVEGHRRGRLLERGTQSGRQS-VKGVLVRKPGAARSARHSPEPLPLQSYPLKRLGGLGGGIKGAPFHPLWTVPCYLGDALSEPAAWCPPSPRSPGSPERKARESAPQRQRNRNRPWRP-EAAGGGGEENNLPKSKAQSFAPPEGEGERRLGASGCVNPAPRRGPSAAAVPTRALAALLPPPAPQPLPTFPLHRPRPHPQYRPRPTLLLPAFLLSLAGWSGELSRGL-HPSINTRGVSNNNRGTSLRTQYGDATAPETPLRYLPRSHPTKGSSCSGTAQSPAPRSSSASRARGHPRAAAAGESGLGAEIRPGFSKSTIPGESDTSSPPFSPACAYPTPRPNTSRDPGFKLYYCF
3'5' Frame 3
SPRREKGVTAT-EGEQSPERRLGDRGGVGSSSNFFQISSLGAGGCSSLCCCCW-KFSSSSQ-KYGCTES-SRS-FLLVKLTLRGIVAGGCWSGAHKAGGSLELKGSWCGNLAQRAVRATVPNLSPCRAIP-SGWVVLVGE-REHPFTPFGQSPAISETHLVNQRLGARRAPAPPGARSAKPGSRPRSGRGIEIGPGALKKPREVAVRKTICQNPRHKVLRHLKEKARGAWALAAA-TPLRAGAPPLRLFPLAP-PLSYPRRHRSPSQPSLSTAPVPTPSTAPVQHSFCQLFFFLSPAGVASSAAGFNTPP-IQGGCQIIIGAPPFALNTEMQLRQRPRCDTSPGATPPRVAAVLEPPRAPLLAVPPHLGREDTRGRPRQGSQGWVRRFALVFPSQRFQENRTHPRLPFPLPAPTQHHVLTPLETQVLNCIIA


How do I create a filtered gene list using expression medians Forgive the simple noob question

I have TPM data of ~50k genes (rows) across ~1k cell lines (columns). In R, I would like to output an "intermediate expression" gene list for each cell line, like:

>head(finalExpressedGeneList, n=3)

Celltype1    Celltype2    CelltypeN
gene1        gene5        gene3
gene3        gene6        gene6
etc          etc          etc

I am defining intermediate expression as "gene expression >= median value within each cell line." I have already replaced all 0 values with NA using:

>data[data == 0] <- NA

I then got medians for each column using:

>colMed <- apply(data,2, FUN = median, na.rm = TRUE)

Now I need to somehow filter the rows from each column and copy the associated gene into a new matrix. I started building the code below:

>if (data[,1]>=colMed [1])
   {
   finalExpressedGeneList <- row.names(data)
   }

But I realized that even if I got it working, I would still have to iterate over all of my columns. I'm sure that there is a cleaner way to do this.

EDIT (adding example data)

#create example data
>j = 5
>M <- matrix(NA,j,j)
>M <- sapply(1:j, function(i) `length<-`(1:(j-i+1), j))
>rownames(M) <- c("gene1", "gene2", "gene3", "gene4", "gene5")
>colnames(M) <- c("cellLine1", "cellLine2", "cellLine3", "cellLine4", "cellLine5")
>view(M)

       cellLine1   cellLine2   cellLine3   cellLine4   cellLine5
gene1          1           1           1           1           1
gene2          2           2           2           2          NA
gene3          3           3           3          NA          NA
gene4          4           4          NA          NA          NA
gene5          5          NA          NA          NA          NA

I would like to get an output of:

cellLine1   cellLine2   cellLine3   cellLine4   cellLine5
    gene3       gene3       gene2       gene2       gene1
    gene4       gene4       gene3
    gene5

because the median of cellLine1 has three genes greater than or equal to the median of that column (3), cellLine2 has two genes that are greater than or equal to the median of that column (2.5), etc...

This below should work. It returns a list,

## example data
j = 5
M <- matrix(NA,j,j)
M <- sapply(1:j, function(i) `length<-`(1:(j-i+1), j))
rownames(M) <- c("gene1", "gene2", "gene3", "gene4", "gene5")
colnames(M) <- c("cellLine1", "cellLine2", "cellLine3", "cellLine4", "cellLine5")
#####
genelist = apply(M,2,function(i)rownames(M)[which(i>=median(i,na.rm=T))])

####
#you can call the genes for cellline 1 using:
genelist$cellLine1

Hope I got it correct!


Proteins that cannot form biofilm? I am trying to build a machine learning training set for bacterial protein sequences that form biofilm, and that cannot. I collected the positive sequences from the GO ontology website but for negative sequences I am not sure which sequences to incorporate into my training set.

Can anyone point me to resources for proteins sequences that are known to be not capable of forming biofilms?

One option is to treat all proteins that are not explicitly known to form a biofilm as proteins that are incapable of doing so, but this will likely result in many false negatives.

You could take a look at one-class classification and positive-unlabeled (PU) learning. Both are techniques specifically designed for these kind of problems where only one of the two classes is known. (In this case: only biofilm forming proteins are known.)


Can I create a CRAM file with a relative reference path? I’m trying to create a CRAM file that stores its path to the FASTA reference as a relative path, rather than an absolute path, so that I can move the files around. Unfortunately I can’t get this to work; I was expecting the following to work:

⟩⟩⟩ samtools view -C -T ../reference/ref.fa -o output.cram input.bam

However, the resulting file contains an absolute path in its header:

⟩⟩⟩ samtools view -H output.cram
…
@SQ     SN:1    LN:249250621    M5:hash     UR:/absolute/path/to/data/mapped/../reference/ref.fa
…

As a result, I am unable to open the file via a different path mount that results in different absolute paths, and I can’t move the file (+ its reference) around, or to different machines.

I know that I could set the REF_PATH environment variable or specify -T when reading the file but I would like to avoid this (the result file needs to be readable by IGV, launched by users who don’t know how to set environment variables).

Is there a way of creating a CRAM file that stores a relative path to its reference?

In your CRAM file, make all of the reference sequences present in the ./ directory. Then make a bash script to make relative links to your reference for whatever mount you are using at the time - ie ln -s TARGET LINK_NAME.

Now your CRAM is portable and you just have to tailor the local working dir to it.


How to convert the .vcf (imputed) file with GT:GP format to GT:DS? I have the genotyped data from impute2 output in .gen format (imputed to 1000G P3). The file has genotype posterior probabilities (GP:3 values per variant). I have converted .gen to .vcf using qctools and the .vcf file has GT:GP format. I need to convert the .vcf file with GT:GP format to GT:DS. Genotype dosages are recommended for use in qtltools/fastqtl analysis. However, I cannot find a tool that would keep the .vcf format and convert GP to DS. Any help much appreciated!

You can do this in Hail.

Here's the rough code to do it (0.1 versions).

Setup:

from hail import *
hc = HailContext()

Import the .gen file. VCF works too:

dataset = hc.import_gen(
    'src/test/resources/example.gen',
    'src/test/resources/example.sample')

Remap the genotype schema and export to VCF:

dataset.annotate_genotypes_expr('g = {GT: g.call(), DS: g.dosage()}')\
    .export_vcf('/tmp/out.vcf.bgz')

Take a look at the getting started page if you want to try it out!

I should note that you may be able to do QTL analyses in Hail, depending on the method you want to use. See blog post here.


How to combine multiple files into one file? I have multiple files (n=86000) with one column each and I want to combine them all into one file with 86000 columns.
I tried the following command
Paste *.txt > mean.txt

And as a for loop
for i in *.txt, do paste $i > mean.txt; done

This did not appear to work because there are 86000 to combine, I then divided my files into multiple folders with each folder not containing more than 15000 files but this did not working
Format of my single file is as follows (It has 16000 lines)
0.047202943506056
0.0383549703366668
0.0159162908608477
0.0238048058550539
0.0694728928685859
0.0389052658410832

What is the bug in the code?

The bug in your loop appears be "> filename" per column, this over-rides the last iteration resulting in a file with a single column, which has been over-written 85999 times. Thus the file output you observe should be the single column present in the last file of your array. If you replaced "> filename" with ">> filename" within the loop you would then append to the file a new column with each iteration. If the loop fell over at some point you could then simply use one of the outputs of wc filename to assess how many columns you had (you might have to transponse the data for wc to work (lots of scripts for this in Perl), although there will almost certainly be an easy solution for counting columns in awk.

Personally, I would use a dataframe for this operation (possibly R, definately Python pandas), because 86000 columns is a lot of columns in any other walk of life, even if it is just post-processing. Generally, if you just want bash, awk is a better solution in this case (I don't code with it, but some here certainly do). However, I do get the idea that paste will work if each column is an identical length.

BTW paste must be lower case, you've got Paste in the first line of your code, but its okay in the second bit of code.

BTW BTW, I did not know shell would take *.txt directly as an array, with Perl/Python you need to "glob" it and assign it as an array/list.


How to specify input resfile for Rosetta? I'm trying to run the following Rosetta command.

rosetta_scripts.linuxgccrelease -s data/binding/dataset/test/protein.pdb \
    -parser:protocol data/binding/test/scripts/pack.xml \
    -out:path:all data/binding/test/output/ \
    -resfile data/binding/test/output/resfile.txt \
    -out:prefix RosettaOutput_ \
    -nstruct 1 -ignore_zero_occupancy false -overwrite

However, I keep getting the error:

core.pack.task.operation.TaskOperations: {0} [ ERROR ] File:resfile.txt not found!

ERROR: Cannot open file resfile.txt
ERROR:: Exit from: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/core/pack/task/operation/TaskOperations.cc line: 1645
Error: {0} [ ERROR ] ERROR: Exception caught by JobDistributor while trying to get pose from job 'RosettaOutput_2MWY_0001'
Error: {0} [ ERROR ] Treating failure as bad input; canceling similar jobs

[ ERROR ]: Caught exception:

File: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/core/pack/task/operation/TaskOperations.cc:1645
[ ERROR ] UtilityExitException
ERROR: Cannot open file resfile.txt

AN INTERNAL ERROR HAS OCCURED. PLEASE SEE THE CONTENTS OF ROSETTA_CRASH.log FOR DETAILS.

protocols.jd2.FileSystemJobDistributor: {0} job failed, reporting bad input; other jobs of same input will be canceled: RosettaOutput_2MWY_0001
protocols.jd2.JobDistributor: {0} no more batches to process...
protocols.jd2.JobDistributor: {0} 1 jobs considered, 1 jobs attempted in 13 seconds

[ ERROR ]: Caught exception:

File: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/protocols/jd2/JobDistributor.cc:327
1 jobs failed; check output for error messages

AN INTERNAL ERROR HAS OCCURED. PLEASE SEE THE CONTENTS OF ROSETTA_CRASH.log FOR DETAILS.

The error resolves itself if I move resfile.txt to the current working directory. Why can't resfile.txt be found?

The problem is my pack.xml specified a resfile.txt location which over-wrote my command-line option.

    <TASKOPERATIONS>
       <ReadResfile name="Design_Resfile" filename="resfile.txt"/>
    </TASKOPERATIONS>

If you remove the filename option in the XML file, you are then able to specify the path.

    <TASKOPERATIONS>
       <ReadResfile name="Design_Resfile"/>
    </TASKOPERATIONS>


How to calculate mutation rate and mutation sites in a genome using FASTA file? I have 6 viral genome sequences of the same virus and 1 reference sequence in FASTA format.

How I can identify mutations and mutation sites in those genomes using FASTA sequences  but how I can do this for FASTA file?
And how I can identify the mutation rate for one genome?

Background If I've FASTQ file then I'll simply align the reads to the reference and by using variant calling tool I will get the mutate sites

You can still align a FASTA file with a tool like bwa mem ( if they are short reads ) or minimap2 for long reads and run it through a variant caller like freebayes. Alternatively, if you have sufficient reads to create an assembly of all of the genomes, or if they already are you can create a MSA of them using a tool like Cactus (as it is built for genomes) or MAFFT if the sequences are short enough and then use your favorite viewer to find variable sites


Using PLINK to find SNPs in LD (Linkage Disequilibrium) with another set of SNPs I have a list of SNPs (index SNPs). As my tagging SNPs, I want to use SNP data from here. The data has the following column headers: Chromosome, Position, MarkerName, Effect_allele,   Non_Effect_allele, Beta, SE, Pvalue.

My ultimate aim is to identify SNPs in the tagging set that are in LD with each of index SNPs. I need to calculate r2 between index SNPs and tagging SNPs within 1MB window and get a table of index vs tagging SNPs (r2>0.8).

How can I use PLINK to carry out this work? I'm new to PLINK and LD clumping

You can't generate LD statistics for summary data from SNPs, as the LD calculation relies on the comparison of variants at the genotype level. Plink expects data to be in a genotype matrix with either each marker for each line (probably better in your case), or each individual for each line. Plink v1.9, in particular, will accept VCF-formatted input files, when populated with 'GT' information that describes the variants for each individual.

If you're okay with using the 1000 genome population for calculating LD, plink v1.9-compatible datasets should be accessible from
here.

If not, you'll need to discuss your project with IGAP to see if you can have access to genotype-level data for LD testing. That's unlikely, because there are personal identification issues associated with non-summary data.

There's comprehensive documentation about the LD command(s) on the Plink 1.9 website here.


full visualisation of draft genomes alignment I have two draft genomes (aprox. 500Mb each), one with ~10'000 scaffolds (genome1.fasta) and one with ~1'000 contigs (genome2.fasta). I would like to visualize their pairwise alignment, which a while ago I found it could easily be done by mummerplot in MUMmer.

In mummer, pairwise draft genomes alignment is supported:

  Many times it is necessary to align two genomes that have not yet been completed, or two genomes with multiple chromosomes. This can make things a little more complicated, since a separate alignment would have to be generated for each possible pairing of the sequences. However, both NUCmer and PROmer automate this process and accept multi-FastA inputs, thus simplifying the process of aligning two sets of contigs, scaffolds or chromosomes.

It is still ok if I do it aligning small portions (e.g. contig vs scaffold), as I could get here. However, if I do it all-vs-all, I end up with a "confetti"-plot (self-coined term) instead of a dot-plot :

I am not sure how to interpret this, and what I would like to see is more of a straight line as the one. Probably it's just a matter of filtering, so I'd like to ask for advice.

I used nucmer through delta-filter as follows:

# mummer version //4.0.0beta1

# nucmer
nucmer --prefix=ref_qry genome1.fasta genome2.fasta

# ref_qry.delta generated from previous step
delta-filter -r -q ref_qry.delta > ref_qry.filter

# plot
mummerplot ref_qry.filter -R genome1.fasta -Q genome2.fasta --filter --layout

Reasons for options -r -q in delta-filter step :

  To determine a one-to-one mapping of each reference and query, combine the options and use delta-filter -r -q.

I'll also accept other methods for visualisation (doesn't need to be withh MUMmer). Any input is greatly appreciated!

EDIT

I found some solutions which I reported as answers. I will accept anyone giving answer on how to tune MUMmer parameters to achieve the goal I describe in the original question.

Found a solution, using D-Genies, worked great.

Some examples from their website:

Thanks to @user172818.


Which tools can detect chimeric RNA (fusion genes) from WGS or RNA-Seq data? Given WGS data or RNA-seq data, which tools can I use to detect gene fusions?

Most of these use RNA-seq data, some use WGS data, and some use both. They are listed alphabetically. I will add to the list when I discover more.

Barnacle: http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-14-550
Bellerophontes: http://bioinformatics.oxfordjournals.org/content/28/16/2114.long
BreakDancer: http://www.nature.com/nmeth/journal/v6/n9/abs/nmeth.1363.html
BreakFusion: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3389765/
BreakPointer: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3561864/
ChimeraScan: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3187648/
Comrad: http://bioinformatics.oxfordjournals.org/content/27/11/1481.long
CRAC: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4053775/
deFuse: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1001138
Dissect: http://bioinformatics.oxfordjournals.org/content/28/12/i179.abstract
EBARDenovo: http://bioinformatics.oxfordjournals.org/content/early/2013/03/01/bioinformatics.btt092
EricScript: http://bioinformatics.oxfordjournals.org/content/28/24/3232
FusionAnalyser: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3439881/
FusionCatcher: http://biorxiv.org/content/early/2014/11/19/011650.full-text.pdf+html
FusionFinder: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3384600/
FusionHunter: http://bioinformatics.oxfordjournals.org/content/27/12/1708.long
FusionMap: http://bioinformatics.oxfordjournals.org/content/27/14/1922
FusionQ: http://www.biomedcentral.com/1471-2105/14/193
FusionSeq: http://www.genomebiology.com/2010/11/10/R104
IDP-fusion: http://nar.oxfordjournals.org/content/early/2015/06/03/nar.gkv562.full
iFUSE: http://bioinformatics.oxfordjournals.org/content/29/13/1700.long
InFusion: https://bitbucket.org/kokonech/infusion/wiki/Home
INTEGRATE: http://www.ncbi.nlm.nih.gov/pubmed/26556708
JAFFA: http://www.genomemedicine.com/content/7/1/43
LifeScope: http://www.thermofisher.com/no/en/home/life-science/sequencing/next-generation-sequencing/solid-next-generation-sequencing/solid-next-generation-sequencing-data-analysis-solutions/lifescope-data-analysis-solid-next-generation-sequencing/lifescope-genomic-analysis-software-solid-next-generation-sequencing.html
MapSplice: http://www.ncbi.nlm.nih.gov/pubmed/20802226
MOJO https://github.com/cband/MOJO
nFuse: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3483554/
Pegasus: http://www.biomedcentral.com/1752-0509/8/97
PRADA: http://www.ncbi.nlm.nih.gov/pubmed/24695405
ShortFuse: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3072550/
SnowShoes-FTD: http://nar.oxfordjournals.org/content/39/15/e100
SOAPFuse: http://www.genomebiology.com/2013/14/2/R12
SOAPFusion: http://www.ncbi.nlm.nih.gov/pubmed/24123671
STAR: http://bioinformatics.oxfordjournals.org/content/29/1/15
STAR-Fusion: https://github.com/STAR-Fusion/STAR-Fusion/wiki
TopHat-Fusion: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3245612/
TRUP: http://www.genomebiology.com/2015/16/1/7
ViralFusionSeq: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3582262/

Other useful programs:

Chimeraviz (visualization tools for gene fusions): https://www.ncbi.nlm.nih.gov/pubmed/28525538 (disclaimer: I created this)

Chimera: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4253834/
OncoFuse: http://bioinformatics.oxfordjournals.org/content/29/20/2539.long
FuMa: http://bioinformatics.oxfordjournals.org/content/early/2015/12/09/bioinformatics.btv721.abstract


How can I use annotations to remove variants not relevant to cancer risk? I currently have ~180 whole germlines and around 10M SNPs/indels. I would like to build a predictive model using Machine Learning (ML) techniques to predict cancer risk according to these germline variants. The thing is, most of these 10M variants are not relevant to cancer, so first I have to annotate them and remove non-relevant SNPs/indels.

Which tools/measurements do you think are important to use in order to filter out irrelevant variants and keep only, let's say, around 10,000 SNPs/indels?

I guess the most obvious is CADD score, so I could keep those SNPs/indels with PHRED >= 10.
Which one would you suggest taking into account that downstream analysis would be ML algorithms that take as features variants that really contribute to cancer?

While your question is specific to cancerous germline mutations, I'd suggest you look at the COSMIC database of somatic mutations to include in your analysis.

There are other factors to include in this kind of analysis you're suggesting, such as predictive deleterious effects (PolyPhen for example can perform such predictions).

If you have 10M variants/indels, for the same form of cancer, then look for common variants, or maybe the frequency of variants identified across exons of a gene.

What you've asked is highly broad and we'd need more specific details about your cohort to suggest more specific approaches.


Multithread fastq processing with kseq.h in C++11? Background

I am using the wonderful kseq.h as my fastq/fasta parser for a C++ project. My actual implementation uses something like this, where I am counting kmers of each sequence and performing a bunch of other processing on every read pair.

Thanks to @user1782818's code my program is very fast, but I want it to be faster. The only way I know how to do this is with multithreading. However, I have never worked with multithreading in C++.

For each read pair I am accessing a vector or map of structs. Within each struct I am incrementing counters or pushing values onto vectors. Essentially I am doing something like this:

//compiled with g++ 7.3 using command g++ this_prog.cpp -o this_prog -lz
// execute with ./this_prog my_file.fastq.gz
#include <iostream>//cout
#include <map>
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char **argv)
{
  //initialize prefix
  std::map<std::string, int> prefix_counter;
  std::string prefix;
  prefix.resize(2);

  gzFile fp;
  kseq_t *seq;
  int n = 0, slen = 0, qlen = 0;
  fp = gzopen(argv[1], "r");
  seq = kseq_init(fp);
  while (kseq_read(seq) >= 0){
    ++n;
    //COUNT PREFIXES
    std::copy(seq->seq.s, seq->seq.s + 2, begin(prefix));
    if (prefix_counter.count(prefix) == 0){
      prefix_counter[prefix] = 1;
    } else prefix_counter[prefix]++;
  }
  std::cout << "Total n. of sequences: " << n << '\n';
  kseq_destroy(seq);
  gzclose(fp);

  //print counts
  for( const auto& sm_pair : prefix_counter ){
    std::cout << sm_pair.first << ": " << sm_pair.second << '\n';
  }
  return 0;
}

Example output from the above program is something like this for a 100,000-read test fastq file.

Total n. of sequences: 100000
AA: 8561
AC: 5178
AG: 7634
AT: 8523
CA: 7526

Question

Q1: Is it possible to perform a task like the one shown above using multithreading and kseq.h to read in a fastq file and process the reads in chunks? I only want to use standard libraries. C++11 is OK - I know that multithreading support was added in C++11.

Q2: Can you please provide a C++ example of how to convert the above or similar program to use multithreading using kseq.h?

You might take a look at my FQFeeder library (https://github.com/rob-p/FQFeeder), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the parsed reads.  The library contains an example program to demonstrate how to use it with multiple threads.  It's worth noting that kseq is already very fast, so multi-threading only really makes sense if worker threads are doing non-trivial work on each read so that you can then make the parsing and processing concurrent (have the processing keep pace with kseq).

I'll note that my library relies on a header-only MPMC implementation and so technically isn't "only using the standard library", but you're already relying on kseq.h, so omitting third-party libraries must not be a hard constraint.


Visualizing protein-ligand interactions with LIGPLOT or any other suitable software? I have used DOCK software to dock 10 ligands into the protein structure.

I want to generate a pdb file that contains all 10 binding modes of the ligands together with the protein. So all together in one pdb file.

Now, what I want to know is when I give this file to LIGPLOT, is it possible to generate something like:

%polar interactions XX
%non-polar interactions XX
%hydrogen interactions XX

Is there a way to do it with LIGPLOT, if now is there some other software that can be recommended for this type of assessment?

I have been having hard time, trying to find a way or software that can perform similar type of calculation and give me similar type of output that I want.

Best

You can visualize with Discovery Studio Visualizer and Schrodinger.
Arpeggio and PLIP can generate output files of the different interactions.


How to interconvert InChI and InChIKey? I would like to retrieve IDs from several databases using InChI as an input, e.g.

InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1

One can use unichem from bioservices for this, however, these functions all require the InChIKey as input, e.g.

KDXKERNSBIXSRK-YFKPBYRVSA-N

Is it possible to interconvert the two using bioservices and if not is it possible to somehow use the functions in unichem with InChI rather than the InChIKey?

I tried:

from bioservices import *
u = UniChem()
u.get_src_compound_ids_from_inchikey('KDXKERNSBIXSRK-YFKPBYRVSA-N')

which works fine, however,

u.get_src_compound_ids_from_inchikey('InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1')

does not work and returns 400.

If it is not passible via bioservices, is there any other python API available for this?

Recently I've developed a python lib for converting InChI to InChIKey: https://github.com/liwt31/chembl_ikey. Maybe it will help.


How can I downsample a BAM file while keeping both reads in pairs? I know how to downsample a BAM file to lower coverage. I know I can randomly select lines in SAM, but this procedure can't guarantee two reads in a pair are always sampled the same time. Is there a way to downsample BAM while keeping pairing information intact?

samtools has a subsampling option:

  -s FLOAT:
  Integer part is used to seed the random number generator [0]. Part after the decimal point sets the fraction of templates/pairs to subsample [no subsampling]

samtools view -bs 42.1 in.bam > subsampled.bam

will subsample 10 percent mapped reads with 42 as the seed for the random number generator.


Highlight CpG sites with Texshade I am using Texshade to produce figures for MSA. Upon reading the documentation I could not find anyway to highlight how conserved CpG sites within the alignment. Does anyone know if this is possible? If not is there another MSA visualizer I can use that can produce a figure like this. Thanks!

The best solution I could find was to use the \shaderegion function on all the sequences so it ended up looking something like this:

\shaderegion{1}{CG}{Red}{OrangeRed}

\shaderegion{2}{CG}{Red}{OrangeRed}

...

Then I cast the rest of shading similarity to grey so the CG pairs were more easily distinguishable.

\shadingcolors{grays}

Two examples one with \fingerprints and one regular.


Alignment for predicting DNA hybridization? I am currently working on a Computer Science project where we are trying to build a large set of orthogonal single-stranded DNA sequences. The goal would be to ensure that when put in solution, the ssDNA do not interact with each other to form double-stranded DNA.

I hope that this makes sense, but my knowledge  biological understanding is developing because I have a CS background. Please excuse any inadvertent errors.

As one of the first steps, we used NUPACK to predict the percent of the resulting AB complex, given two single-stranded sequences A, B in concentration 1μM each. We would use the concentration of AB to classify if the strands are likely to form double-stranded DNA.

However, we cannot use NUPACK for large sets, as we would need to do the above check against each of the other sequences, which is a quadratic problem. For a set of size 30,000 we would have to do 450 million NUPACK computations which is simply unfeasible.

So I turned to the idea of approximating NUPACK. In this paper, it is claimed that

  Ideally, the thermodynamics of hybridization between a candidate
  probe and potential off-target sites would be modeled in
  silico and employed as a means of identifying probe oligos likely
  to only bind their intended targets in a given set of reaction
  conditions. Although powerful utilities such as NUPACK (52–
  54) are capable of performing such simulations, the limited
  throughput of these programs renders a direct thermodynamic
  approach impractical for genome-scale probe design. However,
  we hypothesized that features in rapidly calculated data such as
  alignment scores may be predictive of thermodynamic behavior
  and could therefore serve as a proxy for the information that
  would be produced by thermodynamic simulations.

The alignment scores are then used to predict the interaction. bowtie2 is used in the paper for alignment purposes. However, the context is different as their computation is done on genome information. Very sensitive settings of bowtie2 seem to be slow (hours) on my laptop and result in enormous files (tens of GB). I used less sensitive settings as a compromise but it seems to negatively affect the LDA.

I am looking at alternative applications for alignment, but there are many alternatives (some of the more promising ones so far, in my opinion: SNAP, MMseqs2, Genoogle). Would any of these, or others not mentioned, be better suited for my purpose?

Also any other advice or suggestions are welcome!

As a rough and ready and very easy to perform test I would put the entire universe of dsDNA through Primer3 (details below). This could be done by "splicing" your probe of interest onto each "chunk" of dsDNA you pass through primer3 and restrict the primer search to the 20bp primer you are interested in. In other words you are forcing the program to consider the internal and trans-hybridisation of every possible dsDNA against your probe. If you probe is not recovered from a given search you would flag this on a subsequent parsing script and then that "chunk" of DNA would be assessed further. All the components required are in the link below. You would need to provide your own wrapper, but its fairly easy.

Primer3 will assess the internal thermodynamics of your probe and it has loads of options to assess the full range of behaviour, for example 3' hybridisation (which could be extended because you want 5' to 3'). It would be easy to construct a parameter grid around all the parameters of interest. Ultimately it would provide an additional layer of certainty, and everyone knows and trusts this algorithm (its also written in C, so its fast).

The description of all the components is here , my specific purpose for customising the pipeline is different.

Bowtie2 will take a lot of RAM, I am not sure I would use it for this purpose, particularly if RAM is limiting (its used for genome assemblies against a reference). I'd ask the users whether Blat would be useful, because that is friendlier to the calculation you are performing. Anyway there some heavy duty Bowtie2 users around and their input would be useful.


Does the MAPQ=0 fraction of a BAM file depend on the insert sizes? When doing Illumina 2x150bp sequencing of genomic DNA, and after aligning the reads to GRCh38, does the percentage of the non-N fraction of the human genome as MAPQ=0 depend on the insert sizes of the genomic fragments?

This is, for two identical samples with identical final coverage, sample A having an average insert size of 250bp and sample B having an average of 450bp, would the fraction of MAPQ=0 change between the two?

Yes, as a general rule of thumb mappability increases with insert size (up to a limit) and read length. Whether this will actually occur in a given case will depend more on how randomly the sequencing samples from the genome to begin with (i.e., if the library prep happens to select for/against high mappability regions then the insert size won't much matter). The basic reason behind this is that if one end of a read aligns to a repetitive element that the likelihood of being able to use the other end as an anchor will increase a bit if it's further away. Now having said that, Illumina sequencers in particular like a certain fragment length for efficient sequencing, so you can't go too crazy on increasing the insert sizes and still get decent output (unless you modify your library prep).


How to isolate genes from whole genomes for phylogenetic tree analysis? I have 446 whole Klebsiella Pneumoniae genomes I want to build a phylogenetic tree from. After reading about constructing phylogenetic trees it seems the only option for large numbers of genomes is to isolate a gene with low variability from generation to generation and use this gene to build a tree. For example Lars Jensen recommends using "16S rRNA [or] all ribosomal-protein-coding genes" https://www.biostars.org/p/1930/. What program isolates these genes of interest from the whole genome fasta files and can put them into a multiple alignment file? Or outputs them in a formate ready for a multiple alignment program such as Muave? The reason I say a multiple alignment file is because this it the type of file most phylogenetic tree programs take (I.E. clonalframe).

There are a lot of ways to do this. I suggest using Prokka/Roary to produce a core genome alignment. There's a useful tutorial on the Roary website:

for file in *
do
    prokka --kingdom Bacteria --outdir "${file%%.*}"  --genus Listeria --locustag "${file%%.*}" "$file"
    mv "${file%%.*}"/PROKKA_07052017.gff GFF/"${file%%.*}".gff # use current date
done
roary -f Alignment -e -n -v GFF/*.gff

Alignment/core_gene_alignment.aln can be used as input for phylogenetic analyses


How can I improve ESS values for ucldMean in beast? Using starBeast I ran a 50 locus tree, with 34 individuals, using evolutionary models indicated by AICc in jModel test. I ran this tree with a MCMC chain length of 500 million, storing every 10,000 generations. I used a relaxed clock log normal.

I get low ESS values for posterior, prior, birthRate.t:Species, YuleModel.t:Species, popMean, and TreeHeight.Species.
ucldMean.c:Locus_name is low for each loci.
The TreeHeight statisitcs also have low ESS values.
Also rate.c:Locus_name.mean for each loci is low as well.

Is there a way to change my parameters in Beauti to help increase the ESS values of these statistics? Is there a parameter that would likely cause the ESS values of these particular statistics to be low? Would increasing the MCMC chain length increase the ESS values, given that I've already run half a billion?

The official BEAST website suggests a number of solutions to improve the ESS of a parameter

It is not uncommon to have several billions of generations so I don't imagine increasing your number of generations would hurt.


single stranded cDNA vs double-stranded cDNA Single stranded cDNA is used by race pcr and rt-pcr while next-gen sequencing (Illumina) requires double stranded cDNA. What is the reason behind this ?

Sequencing also requires single stranded DNA, otherwise the sequencing primers can’t bind. Sure there are steps in the protocol where double stranded material exists, but that’s transient.


Is there a tool that can take a protein's amino acid sequence and would display it's locus on the genome? I have the UNIPROT IDs, PDB IDs and FASTA files of several known proteins. I am looking for a tool that can take as input the protein's amino acid sequence and display the coding nucleotides of those amino acids. Moreover, I look for a tool that would be able to display those nucleotides' location on the human genome.

The final outcome I desire to achieve is to find SNPs in the exons of the protein.

A tool that can handle several amino acid sequences at once is preferred for I have plenty of proteins. I've tried to use UCSC's BLAT human genome search but it was too cumbersome. I've also tried BLAST, but I'd like a more "straightforward" method for displaying the genome locus for each protein. By the way, I've tried to use UCSC's Table Browser but couldn't query several genes simultaneously.

Thank you for your assistance.

I think the mistake here is thinking that you need to use the sequence in order to find the genome co-ordinates. If these proteins come from a well annotated speices, you can easily do this with BioMart.

Go the the ENSEMBL web site.
From the top menu, select BioMart
Under "-CHOSE DATABASE-" select Ensembl Genes, and under "-CHOOSE DATASET-" select your species (eg. Human Genes).
Select "Filters" from the lefthand menu and then expand the "GENE" section.
Select the "Input external references ID list [Max 500 advised]" and from the drop-down select "PDB IDs" or the correct "Uniprot KB" ID.
Paste your IDs into the box.
Select "Attributes" on the left-hand menu. Making sure you've the "Fetaures" radio box select, you can select Chromosome and either Gene Start/Gene End or Transcript start/Transcript end from the GENE section. These coordinates will also contain the untranslated regions unfortunately.  You may also want to play about with the various things under "PROTEIN DOMAINS AND FAMILIES", particurlarly under "Protein Features", although I'm not sure what they all mean.
When you are happy with the selection, select the big "Results" at the top to get the coords.

You can also get out Sequences for both the full transcript and the CDS and known germline and somatic variants that fall within the boundaries of the genes.


Using column 2 of one file to match with two columns of another file, and append I have file 1 like following:

1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    16855676    16855696    CTGTCCAGCAGAGGGCGGTG

And file 2 as following

1   15776231    2   5008    G:5002  A:6 1   16855677    2   5008    A:5003  C:5

I am trying to write a code(but failing!) to do the following:

if $2 of file2 comes in between $2 and $3 of file1 (between the interval of $2 and $3) and if $8 of file2 comes in between $5 and $6 of file1 (between the interval of $5 and $6), then the output would be $1,$2,$3,$4 from file1 and $2,$3,$4,$5,$6 from file2 and $5,$6,$7 from file1 and $8,$9,$10,$11,$12 from file 2, all in one line. So something like this:

1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776231    2   5008    G:5002  A:6 16855676    16855696    CTGTCCAGCAGAGGGCGGTG    16855677    2   5008    A:5003  C:5

Something like bedtools does but I am not able to use bedtools one this since I don't have two columns in file2. Is there any way to make bedtools work on this?

(There could be many lines from file2 intersecting between the same line in file1,so that line could be repeated)

Pipe a modified form of the second file to BEDOPS bedmap and the first file, then pipe that result to cut out desired columns:

$ awk -v OFS='\t' '{ print $1, $2, ($2+1), $3, $4, $5, $6, $7, $8, $9, $10; }' second.txt | bedmap --echo --echo-map first.bed - | cut -f1-4,8,9,11- - > answer.bed

This should run pretty quickly and use very little memory, as it uses Unix streams on (presumably) sorted BED data.

You'll want to test the column indices in my answer by running this on your test inputs, but that's an easy adjustment for the larger dataset, in case I missed a column.

If you're not sure about sort order, pipe to sort-bed:

$ sort-bed first.unsorted.bed > first.bed
$ awk -v OFS='\t' '{ print $1, $2, ($2+1), $3, $4, $5, $6, $7, $8, $9, $10; }' second.txt | sort-bed - | bedmap --echo --echo-map first.bed - | cut -f1-4,8,9,11- - > answer.bed

Use streams where you can. It's a different way of thinking but can pay huge dividends.


how to calculate coverage of each base of mapped reads? I have a RNA aligned file without any header info.
I need to calculate the coverage of each mapped base for per chromosome. Is it possible to do it with perl ?

NS500377:28:H16A7BGXX:2:21205:16038:7658        99      chr2   3000482 255     46M     =       3000502 56      GTCCTACAGTGTGCATTTCTCGTTTTCCACGTTTCT    AAAA)FFFFFFFFFAFFFFFFFFFFF.F.FFFFFAF    XA:i:2  MD:Z:21T4T7T1   NM:i:3
NS500377:28:H16A7BGXX:3:11503:24953:10427       99      chr2   3000482 255     50M     =       3000513 71      GTCCTACAGTGTGCATTTCTCATTTTTCAAGTTTTTCAGT        AAAAAFFFFFFFFFFFFAFFFFFF7AAFFFFFFFFFFFFF
        XA:i:1  MD:Z:21T7C10    NM:i:2
NS500377:28:H16A7BGXX:3:12404:23929:9392        163     chr2   3000482 255     49M     =       3000483 40      GTCCTACAGTGTGCATTTCCCTTTTTTCACGTTTTTTAG AAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:1  MD:Z:19T16C2    NM:i:2
NS500377:28:H16A7BGXX:3:12404:23929:9392        83      chr2   3000483 255     49M     =       3000482 -40     TCCTACAGTGTGCATTTCCCTTTTTTCACGTTTTTTAGT AFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:2  MD:Z:18T16C3    NM:i:2
NS500377:28:H16A7BGXX:3:21606:15675:1639        83      chr2   3000483 255     50M     =       3000480 -43     TCTTACAGTGTGCATTTCTCTTTTTTCACGTTTTTCAGTG        AFFFFFFFFFAFFFFFFFFFFFFFFFFFFFFFFFFAAAAA
        XA:i:0  MD:Z:2C37       NM:i:1
NS500377:28:H16A7BGXX:4:12605:18505:11088       83      chr2   3000486 255     50M     =       3000481 -45     TACAATGTGCATTTCTCTTTTTTCACGTTTTTCAGTATTT        FFFFFFFAFFF)FFFFFFFFF<FFFF<AFFFFFFFAA<AA
        XA:i:1  MD:Z:4G31G3     NM:i:2
NS500377:28:H16A7BGXX:2:13303:20634:12085       99      chr2   3000487 255     49M     =       3000521 73      ACAGTGTACATTTCTCATTTTTCACGTTTTTCAGTGTTA AAAAAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF. XA:i:2  MD:Z:7G8T21T0   NM:i:3

You are looking for bedtools genomecov, specifically -d option. See https://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html , for example:

bedtools genomecov -i A.bed -g my.genome -d

I recommend not reinventing the wheel: use bedtools instead of reimplementing the same functionality in Perl.


Making PLINK compatible files from VCF file without phenotype information I have a big VCF file that I need to convert to, preferably, bed/bim/fam files that are readable by plink. Currently using plink 2. I am aware that this version of plink can be used to convert VCFs into binary peds using
plink --vcf file.vcf.gz --make-bed --out out
Or something similar. I also know that in this case, plink will automatically fill the phenotype and sex columns in the bed and fam files with zeros. I know there is a --pheno flag that allows specifying a txt file for the phenotypes and a --make-pheno flag that can be used for the same effect.
How exactly does this work? How can I encode the phenotype and sex information correctly? Can I just use --make-bed and then manually make a phenotype/sex file to use down the line with the --pheno flag? Or would I do something like
plink --vcf file --pheno phenotype_file --make_bed --out out_file
directly? In this case, will there not be a problem with the fact that the generated .bed in the first step does not have this information? What about with sex?
Another option I can think of is just edit the generated .fam to include the correct information and use it downstream, but that information will still not be in the .bed file.
I have no family relationships, so the FID fields would just be 0s (I guess I can add the --const-fid flag for that effect when I'm converting, right?).

You can use --pheno to specify a phenotype file you want to use, the 1st and 2nd columns being FID and IID. These are just used on the fly and cannot be written to the BED files as I understand. --make-pheno is used to write binary phenotypes, which is unsuitable if your phenotype is actually continuous.
But I still do not quite understand what you want to do here.


Can I stop my nanopore sequencing run if there are no more reads being produced? I am sequencing a whole genome on MinION. I have used this flow cell for 6 hours for the phage lambda control. That gave me around 300,000 reads. I started a 48h sequencing run on the same flow cell after washing it. The sequencing has been going on for 1 day and 21 hours. It is not producing any reads at the moment (Strand = 0). The number of single pores is also 0. Can I "Stop acquisition" without negatively affecting the results?

The image below is a screenshot taken from MinKNOW, the software used to control the MinION.

Yes, you can click "stop acquisition" and your run won't be negatively affected. All of the reads are saved as they are generated. I am not sure how this will impact live basecalling though if that is something that you do.


read and matching pattern with python I have read the content of a text file into pandas and needed some help matching the pattern.

Here is the pattern, where a can any number greater than zero and b can be zero or any number greater than zero:

father=a|0 || 0|a
mother=a|0 || 0|a
daughter1=b|0 || 0|b
daughter2=a|0 || 0|a
daughter3=b|0 || 0|b
son1=a|0 || 0|a
son2=b|0 || 0|b

Here is the function that reads the content of text file to pandas

import pandas as pd
def read():
    varian=pd.read_csv('file.txt',sep = "\t")

father=varian["FATHER"]
mother=varian["MOTHER"]
daughter1=varian["DAUGHTER1"]
daughter2=varian["DAUGHTER2"]
daughter3=varian["DAUGHTER3"]
son1=varian["SON1"]
son2=varian["SON2"]

The content of the file looks like this

HROM    POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3   SON1    SON2    INFO    FREQUECY
1   1226852 G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=12   0.005494505
1   1847936 C   T   0/1 1/1 1|1 1|1 1|1 1|1 0|1 AC=1;AN=2184    0.000457875
1   2428427 C   G   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AC=4;AN=2184    0.001831502
1   2515616 G   A   1/0 1/0 1/0 1|1 1|1 1/0 1/0 AC=1;AN=2184    0.000457875
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755

This is not a bioinformatics question, but a basic pandas question, so SO is a better place for this... But your odd looking VCF file can be parsed as follows without even using regex:

import pandas as pd

def splitter(cell, key):
    parts = cell.split('||')
    if key in parts[0]:
         p = parts[0].split('|')
    else:
         p = parts[1].split('|')
    sp = [x.strip() for x in p if key not in x]
    return float(sp[0])

def read():
    splitter_a = lambda v: v.split('||')[0].split('|')[1])
    varian=pd.read_csv('file.txt',sep = "\t")
    a_splitter = lambda c: splitter(c, 'a')
    b_splitter = lambda c: splitter(c, 'b')
    varian = varian.assign(father_a=varian.FATHER.apply(a_splitter)\
              .assign(father_a=varian.FATHER.apply(b_splitter)\
              .assign(mother_a=varian.MOTHER.apply(la_splitter)\

etc.

display(varian) or print(varian.columns) will show that you have new columns called father_a etc.


Tools or softwares for clustering evaluation I am using a few clustering algorithms as well as my own tool for protein sequences. Also, I have benchmarked clustering results (i.e. ground truth) to compare with the results. Is there any tool/software/script to generate evaluation metrics for the clustering results.

The clustering results/benchmarked are in the following format

1  Name1
1  Name2
2  Name3
2  Name4
2  Name5
.....

UPDATE: I am evaluating the clustering performance using scikit

The definitive answer is no, that is why they are not so popular. There are metrics that can be used:

you can examine "minimum distance" as a criterion under the assumption that the lowest distance between taxa represents the largest homology (e.g. if it is sequence data), which would be translate to other biological scenarios.
resampling, e.g. bootstrapping, jack-knifing and examine the tree structure of each node against the biological plausibility of the answer being correct.

For both points a biological a priori, i.e. biological plausibility of assessing the tree structure is required to confirm/reject the clustering method.

Clustering algorithms are not like the probabilistic methods, where the highest likelihood is the singular criteria and of course they do give different answers. Its a brief explanation I know but hope it makes sense.


Total reads aligning to each reference within a bam file I have two PCR amplicons that have been multiplexed and sequenced using the nanopore minion.

I have aligned the fastq reads using minimap2 with a reference file containing both amplicon sequences and generated a bam file that I have viewed using IGV.

I am looking for a way to generate some simple summary statistics.

In particular, is there a way to extract the total number of fastq reads aligning to each amplicon reference from the bam file?

The quick way to get the number of alignments on each reference is

samtools idxstats my_bam.bam

Number of reads on each reference is column 3. Although, as has been pointed out, this will give you the total number of alignments per reference, not the total number of reads (each read might give rise to more than one alignment). That said I do tend to us this as generally I'm after a rough approximation, rather than an accurate number.

In theory, only one alignment for each read should be marked as primary, so the following should give you what you need quickly and at low memory usage:

samtools view -bF 2304 my_bam.bam > primary_only.bam
samtools index primary_only.bam
samtools idxstats primary_only.bam


What is the purpose of 'Folder locations' in MaxQuant? MaxQuant is a tool for proteomics quantification. In the settings under Global parameters / Folder locations are three folder locations that are possible to be changed. Thought the manual does not even mention these options. Is there any benefit in changing these locations? Or can it be bad to do that? For example, if I run multiple MaxQuant instances at the same time doing different searches?

What is the meaning of the Fixed andromeda index folder, Temporary folder and Combined folder location in MaxQuant? What is it used for and when would you (not) change it?

These folders are usually generated in a default directory (where the raw files are). It can however be beneficial to set the paths to these folders manually.

  Temporary folder

that is obviously where temporary files go, so having an SSD or other fast drive will benefit the I/O operations. Set that path to a fast storage partition.

  Combined folder

all the output, including txt directory will go here. Setting it to fast storage might speed up writing the txt files. If you set it to the same path while running multiple instances of MQ it will obviously override output.

  Andromeda index

this is where index files from Andromeda search are. By default these are inside a combined/search folder. Since these index files can get big, once generated by MQ they can be reused. Set this path if you want to reuse previously generated index files.


Why indexing sequence file with minimap2? I have been using minimap2 for a few month now and I just found a code that follow these steps :

minimap2 [-x preset] -d target.mmi target.fa
minimap2 -a [-x preset] target.mmi query.fa > output.sam

when I never saw this .mmi index. I always went directly like this :

minimap2 -a [-x preset] target.fa query.fastq > output.sam

So my question is why should I index my reference genome before mapping ? Meaning is it only for speed up ? Will it up speed my code if I do it ? Or it is already done without specifying it ?

If you don't provide an index, minimap2 will just generate the index on the fly (so, yes, it's done without specifying it). As such the mapping step will take longer. However, index generation is supposed to be very fast in minimap2 (a few minutes for the human genome), so it probably won't matter much.

Doing it separately will therefore also not speed up your code unless you end up using the pre-generated index multiple times.


What to do with the configuration test cell after configuring the MinION with it? What to do with the configuration test cell after configuring the MinION with it?

I received a MinION with a configuration test cell in it. After running a configuration in the MinKNOW software while the CTC was in the MinION, I wonder what to do with the CTC? Am I going to need it later or can I throw it away? If I should keep it, in what conditions should I keep it (temperature, light, etc.)?

Also, do you know why we need to do a configuration of the MinION the first time?

I posted this question on Biology Stackexchange as well but maybe the community here will be more receptive?

I would keep it, in case you need to do another configuration later on. I have used mine for that purpose.

I think storage doesn't really matter, I just put it in a drawer somewhere. Keep it away from liquids and extreme temperatures :)


How do I find identical sequences in a FASTA file? I want to create a database for a proteomics study. Therefore, the mapping from a given sequence to a protein ID has to be unique. I am wondering whether there is already a built-in function in Biopython for that, but I could not find any. The duplicated sequences should be merged into a single entry and there respective IDs should be just concatenated (see expected outcome below).

Using the following example input

>prot_id_1
MAWIGLDISKLFVENSRDAAA
>prot_id_2
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_5
MAWIGLDISKLFVENSRDAAA
>prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_7
MAWIGLDISKLFVENSRDAAA

I currently do it as follows:

from Bio import SeqIO
import numpy as np
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC

# Read entries into dictionary, gi: sequence
record_dict = SeqIO.to_dict(SeqIO.parse("test.fasta", "fasta"))

# Extract sequence as string
record_dict = {gi: str(record_dict[gi].seq) for gi in record_dict}

# Check for duplicated sequences
seq_id_dict = {}
for k, v in record_dict.items():
    seq_id_dict.setdefault(v, []).append(k)

# Filter sequences that appear more than once
filter_dict = {seqi: li for seqi, li in seq_id_dict.items() if len(li) > 1}

# Merge identifiers to a single one
merged_ids = {"_".join(sorted(li)): seqi for seqi, li in filter_dict.items()}

# Remove all entries that are associated with duplicate IDs
for gi in set(np.concatenate(filter_dict.values())):
    record_dict.pop(gi, None)

# Add the merged IDs again
record_dict.update(merged_ids)

# Create list of Seq objects
final_seq = [SeqRecord(Seq(seqi, IUPAC.protein), id=gi, name='', description='') for gi, seqi in record_dict.items()]

# Write file
SeqIO.write(final_seq, 'test_output.fasta', 'fasta')

which gives me the desired output:

>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS

Is there a smarter way of accomplishing this?

The trick would be to swap the key in the dictionary to be the sequence itself.

Also I would recommend using a different separator that "_" since that is what the current ids have so that you can easily separate the individual ids from the concatenated id. I used a pipe "|" in this example. Also I just manually wrote the FASTA output instead of using Biopython to do it.

from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")


How to predict protein binding to ligands? I want to predict from protein sequence if the protein binds to metal, nuclear or small ligand. How can I do this ?  Which features are relevant if I want to use them in a machine learning algorithm ?

If you are training a machine learning algorithm you probably want to train against things you know to be true, rather than making predictions from the sequence and training against the predictions.

In this case you will want to use protein structures with the relevant ligands bound as your data. If you go to the RCSB PDB advanced search you can search for entries with any ligand present, specific ligands present, protein and DNA present, etc. Then you can extract the FASTA sequences for these entries.

See also databases like MetalPDB for specific ligand types, and one of our papers for an example of creating a dataset of metal binding proteins.

Features to use would be the raw sequence as a one-hot encoding on amino acid types. You could add other sequence level features such as secondary structure predictions, disorder predictions or physico-chemical properties. However a decent network would be able to predict much of this if it needs to. You could create alignments and put in a sequence profile or covariation data, which would be giving the network more information.


CRISPR Sequence Finder and Database Download I am searching for tools to pull CRISPR Spacers from Bacterial Genomes. I am aware of the CRISPRDB and the corresponding identification tool on the web server.

Are there other tools for finding CRISPR Spacers?
Is there a way to download the spacer sequences from CRISPRDB?
Is there another database that contains all validated / inferred bacterial CRISPR sequences?

CRT (CRISPR Recognition Tool) http://www.room220.com/crt/

Crisprs Finder Online Tool https://crispr.i2bc.paris-saclay.fr/Server/

Piler-CR http://www.drive5.com/pilercr/

(Omic Tools List) https://omictools.com/crispr-detection-category


Nonbonded interactions on the GPU i'm using this command
gmx mdrun -deffnm nvt0 -nb gpu

but got
Nonbonded interactions on the GPU were requested with -nb gpu, but the GROMACS
binary has been built without GPU support. Either run without selecting GPU
options or recompile GROMACS with GPU support enabled

I've Nvidia rtx2060 with cuda 11.1 on ubnutu 18.4.5

As the error tells you: the GROMACS binary has been built without GPU support.
Did you install gromacs with apt-get? Usually these kind of programs are shipped with a basic built, with a lot of features disabled.
If you did built gromacs from source, then cuda was not found during configuration.
I would try building gromacs from source and:

be sure cuda is in your path and check that nvcc and the drivers have a matching version
run cmake and check the output, to see if the video card is recognised
run cmake with -DGMX_GPU=ON

Note that it seems that gromacs doesn't support multiple gpu accelerations( e.g. it is either openCL or cuda, you can't have both)


How to create Phylogenetic Trees from fasta files in Python or R? I have around a hundred Fasta files (and will collect several thousand) with DNA sequences and +50x coverage. What is a recommended method to construct a phylogenetic tree? Solutions in Python or R are sought.
I found Phylo from Biopython only handles already calculated trees.

I would not look for a package for this, but instead build a small pipeline calling external tools with something like the following workflow:

Cluster the ~100 sequences with CD-HIT-EST/PSI-CD-HIT or many other options
Take all the sequences that form one individual cluster and build a multiple sequence alignment (MSA) with MAFFT/ClustalOmega or similar
Take the MSA and build a phylogenetic tree with a Maximum-Likelihood approach like iqtree or similar
Visualize the tree file with Jalview or similar

Of course this is rather general and depending on exactly what you're doing you may want a different workflow and/or different tools. You should also explore the parameter space, do not assume the defaults are necessarily good choices


How is BLAST's nr database created? Is there a paper or web page describing the procedure for creating the nr database used by NCBI's BLAST implementation?

I presume it's some type of clustering, but I'm curious about how exactly sequences are condensed into non-redundant representatives.

The Refseq team and also the NCBI resource coordinators team publish a new paper every few years, so check out the many papers (e.g. here or here), but to answer your 2nd question, non-redundancy here is (I think) defined very strictly as proteins that are identical in terms of sequence and length, so the clustering is trivial, without the need for a sophisticated clustering algorithm as required to detect more remote homologs.


How to extract unmatched reads using bwa and samtools? I have a single read (NOT paired) that I need to pass through the workflow described in Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.
Here in Materials and methods the procedure is described as the following:
The first step of the workflow consists of an alignment of the reads against the host genome (Fig. 2i).
This step aims at discarding reads that map to the host genome and may partially map to the viral genome after segmentation, and at reducing the working file size.
For example, MV and rMV-ΔV data sets were generated from total RNA samples of infected cells and contained mostly reads mapping the human genome (≈99%).
This step uses a combination of bwa mem and samtools view with the parameters –bS –f4.
An additional step consists of an alignment of the reads against the viral genome of interest, in order to exclude reads perfectly mapping the viral genome.
Therefore, only unmapped reads are further analyzed.
Of note, clipped reads (i.e., CIGAR motif contain S or H) are also conserved. Some of these reads may map to viral genome recombination junctions that are present in DI genomes.
I have already been suggested to use bowtie2 instead of bwa but, firstly, the output is not clear for me and secondly, I would like to test the official protocol.
Since the article suggests to use bwa and samtools on this very first step, that is what I have done so far:

(Optional) not sure if this is important but as suggested by someone I have transformed .fna --> .fa
cp GCF_000001405.39_GRCh38.p13_genomic.fna GCF_000001405.39_GRCh38.p13_genomic.fa

Indexed human genome with bwa
bwa index GCF_000001405.39_GRCh38.p13_genomic.fa

aligned my only read with indexed genome from step 1
bwa mem GCF_000001405.39_GRCh38.p13_genomic.fa /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz > whole.sam

Transform .sam --> .bam
samtools view -S -b whole.sam -o whole.bam

Separated unmapped reads (as it is recommended in Materials and Methods using -f4)
samtools view -f4 whole.bam > sample.unmapped.sam

Converted unmapped reads into .fastq format (since this is the format used by the software later)
samtools fastq sample.unmapped.sam > unmatched.fastq

I wonder if these steps are correct?

Yes, correct but overly many steps. There is no need to convert between bam and sam. samtools can read from stdin and handles both sam and bam and samtools fastq can interpret flags, therefore one can shorten this to:
bwa mem (...options) | samtools view -o out.bam
samtools fastq -f 4 out.bam > unmatched.fastq


Can I pilon-polish long reads with Illumina short reads to improve structural variant detection? I have pacbio Sequel data at 50 x coverage for a strain of animal. I would like to find structural variants compared to the reference genome sequence. At the moment, I align my reads to the reference genome using NGM-LR and call the variants using Sniffles. The resulting VCF file contains variants that are well supported in a visual inspection in IGV as well as some more questionable variants (especially larger events). I have used the reads corrected by canu for alignment. I think it splits my long reads into very long reads and not-so-long reads and uses the shorter reads to correct the longer reads. However, in IGV, my reads still have a considerable number of mistakes compared to the reference. I was wondering if it is OK to use pilon to polish my long reads with some Illumina data from the same sample? The main purpose of pilon is to correct assemblies but I was wondering if it can be used to correct subreads from PacBio sequencing?

I'd recommend using Canu for correction rather than Pilon, because it has a component that is specifically designed for read-level correction. The newest version of Canu (v1.7) will track all reads through the correction and assembly process:

  Ensure that every raw read is either corrected or used as evidence for correcting some other raw read. This serves to rescue short plasmids in high coverage datasets, and it is no longer necessary to set corOutCoverage to achieve the same result.

More details here.


How can I plot multiple sequence alignment in R as a heatmap I have many amino acids sequences in a fasta format that I did multiple sequence alignment. I was trying to plot something like binary a code as heatmap. If it had a change in the AA compared to the reference sequence it would be red, if it did not change would be yellow, for example.

I came across to msaplot from ggtree package. I also checked ggmsa package for that. But so far, I did not get what I wanted. So basically I wanted to:

change the multiple sequence alignment to a binary matrix (if the amino differs from the reference sequence, plot x, if not y)
plot as a heatmap with two colors based if the amino acid changed or not

here's some data example of what I doing, just instal ggmsa and you can try this:

protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
ggmsa(protein_sequences, start = 265, end = 300, color = "Chemistry_AA",font = NULL)

@StupidWolf answered me on stackoverflow and he's answer solved my question. I will share it here as well.

@StupidWolf:

We read in the alignment:

library(Biostrings)
library(ggmsa)
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
aln = readAAMultipleAlignment(protein_sequences)
ggmsa(protein_sequences, start = 265, end = 300)

Set the reference as the 1st sequence, some Rattus, you can also use the consensus with consensusString() :

aln = unmasked(aln)
names(aln)[1]
[1] "PH4H_Rattus_norvegicus"

ref = aln[1]

Here we iterate through the sequence and make the binary for where the sequences are the same as the reference:

bm = sapply(1:length(aln),function(i){
as.numeric(as.matrix(aln[i])==as.matrix(ref))
})

bm = t(bm)
rownames(bm) = names(aln)

The plot you see above has sequences reversed, so to visualize the same thing we reverse it, and also subset on 265 - 300:

library(pheatmap)
pheatmap(bm[nrow(bm):1,265:300],cluster_rows=FALSE,cluster_cols=FALSE)

The last row, is Rattus, the reference, and anything similar to that is read, as you can see in the alignment above, last 4 sequences are identical.


Publicly available genome sequence database for viruses? As a small introductory project, I want to compare genome sequences of  different strains of influenza virus.
What are the publicly available databases of influenza virus gene/genome sequences?

There area few different influenza virus database resources:

The Influenza Research Database (IRD) (a.k.a FluDB - based upon URL)

A NIAID Bioinformatics Resource Center or BRC which highly curates the data brought in and integrates it with numerous other relevant data types

The NCBI Influenza Virus Resource

A sub-project of the NCBI with data curated over and above the GenBank data that is part of the NCBI

The GISAID EpiFlu Database

A database of sequences from the Global Initiative on Sharing All Influenza Data. Has unique data from many countries but requires user agree to a data sharing policy.

The OpenFluDB

Former GISAID database that contains some sequence data that GenBank does not have.

For those who also may be interested in other virus databases, there are:

Virus Pathogen Resource (VIPR)

A companion portal to the IRD, which hosts curated and integrated data for most other NIAID A-C virus pathogens including (but not limited to) Ebola, Zika, Dengue, Enterovirus, and Hepatitis C

LANL HIV database

Los Alamos National Laboratory HIV database with HIV data and many useful tools for all virus bioinformatics

PaVE: Papilloma virus genome database (from quintik comment)

NIAID developed and maintained Papilloma virus bioinformatics portal

Disclaimer: I used to work for the IRD / VIPR and currently work for NIAID.


Improve scRNA-seq dataset for further analysis I got a dataset from C.Elegans scRNA-seq paper:

GSM2599701_Gene.count.matrix.celegans.cell.Rdata in GSE98561_RAW.tar

The dataset is 40 000 x 68 000, where rows represent genes and columns - cells. So, I took it and tried to process myself to build an scRNA-seq pipeline. Here is what I did:

I filtered out the genes that have zero counts in all of the cells and the dataset was reduced to 29 000 x 68 000
I removed all of the cells with counts < 100 in all of the genes - the dataset became 29 000 x 66 000
Then, because the dataset was too big to run normalization even on the cluster with 120Gb RAM (because there are multiple distinct types of cells, first clustering needs to be done), I selected just even columns and ran normalization with 29 000 x 33 000 dataset (UMI_count):

library(scran)
library(scater)
sce <- newSCESet(countData=UMI_count)
clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, clusters=clusters, positive=TRUE)

After running the code above I decided to check whether the data is fine and so I ran:

> summary(sizeFactors(sce))
Min.    1st Qu.  Median  Mean   3rd Qu.    Max.
0.0000  0.0000  0.0000  0.0717  0.0000 33.3900

I also ran PCA on the normalized dataset and it looks like that:

It seems to me that the normalized dataset is terrible and I need to do some further processing before I could run further analysis. What else could I do to improve it? How to filter it? There are no spike-ins and maybe 200 mitochondrial genes. The approach described here does not work, probably because the majority of the cells have low number of genes expressed:

I tried removing low-abundance genes after normalization, but it seems that most of them are going to be removed:

>ave.counts <- rowMeans(counts(sce))
>keep <- ave.counts >= 1
>sum(keep)
109

Should I filter all of the cells - columns - for up to 500 total gene expression count instead of 100? Is it a good idea? I can't think of anything else.

I don't think you can conclude that the dataset is terrible based on that PCA. Depending on the specific protocol, each scRNA-seq dataset is going to be very different. Unlike bulk RNA-seq where all the samples are going to be of very similar quality, individual cells will be highly variable. For example, the most basic QC metric is the number of reads. If some of your bulk RNA-seq samples have 10X more reads than other samples, that's a cause for concern. In scRNA-seq experiments, that's not at all surprising. The PCA exposes the major sources of variation. In scRNA-seq experiment, those are probably going to be technical rather than biological, so the output may not make much sense. Furthermore, even if all the cells are of perfect quality, PCA is not really optimal for sparse data with many dropouts like you see with scRNA-seq.

Regardless of PCA, those library sizes seem very small. Based on your plot, there are barely any cells with over 10K reads. I understand this is C elegans, but that seems somewhat low. That may be why you are having trouble with normalization. It would help with understanding the data if you changed the X-axis on those plots so that you can actually see where the peak is. Right now, it looks like 0, but that can't be the case.


What is the name of this type of figure? This figure comes from this wiki page.

I googled "gene location figure" and "gene location plot", none gets similar results. I also tried search by image, none of the results is close.

so, what is the name of this type of figure?

This kind of visualization showing the banding patterns of chromosomes are called "chromosome ideograms".

You can use for example IdeoViz in R to generate them.


WGCNA correcting p value for varying sample numbers in each of the trait columns I am correcting a data frame of sample comparison using p value. However in some columns I have different sample numbers. Within the data frame I have 4 comparisons e.g B_R, R_C, B_C, C_C. They each have different sample numbers in the comparisons. I need to create a loop or change the code so I can vary the nSamples numbers for each of the comparisons. I have tried to do this for each comparison separately but to no avail. Does anyone have any suggestions? Thank you. I have tried..

nGenes = ncol(cells)
nSamples = nrow(cells)

MEs0 = moduleEigengenes(cells, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor= cor(MEs, status, use="p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")  ``````````

I have tried to code for the varying Sample numbers using:

MEs0 = moduleEigengenes(cell, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor= cor(MEs, NR_status, use="p")
moduleTraitPvalue$B_R = corPvalueStudent(moduleTraitCor, 32)
moduleTraitPvalue$R_C = corPvalueStudent(moduleTraitCor, 31)
moduleTraitPvalue$B_C = corPvalueStudent(moduleTraitCor, 47)
moduleTraitPvalue$C_C = corPvalueStudent(moduleTraitCor, 31)

sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue$NR_R,
                      moduleTraitPvalue$R_C,
                          moduleTraitPvalue$C_C,
                      moduleTraitPvalue$NR_C, 1), ")", sep = "")

You can use the function corAndPvalue. Simply replace the lines

moduleTraitCor= cor(MEs, status, use="p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

by something like

cp = corAndPvalue(MEs, status)
moduleTraitCor = cp$cor
moduleTraitPvalue = cp$p

corAndPvalue takes automatically into account the number of valid samples for each element of the resulting correlation matrix.


How to calculate Fst from AMOVA I calculated an AMOVA from a genind object, with one hierarchical factor. In the table I obtain there are SSD values (for my grouping factor,"Error" and total) and sigma2 values (for my grouping factor and "Error")

I have two questions:

What does "error" stand for?
How do I calculate Fst? Which values do I have to use?

A genind object in R is an object which contains allelic information about a set of individuals.
This is the table that I get:

Analysis of Molecular Variance

Call: amova(formula = plantula.dist ~ g, nperm = 100)

           SSD         MSD  df
g     354992.9 354992.9297   1
Error 310814.9    693.7834 448
Total 665807.9   1482.8683 449

Variance components:
       sigma2 P.value
g     1574.69       0
Error  693.78

Variance coefficients:
       a
224.9956

Error is the estimation that can't be fitted by your parameters (in this case g).
The amova function in pegas package does the following:

  This function performs a hierarchical analysis of molecular variance as described in Excoffier et al.
  (1992). This implementation accepts any number of hierarchical levels

And as per the definition of Fixation index, you need the variance in the frequency of the allele between different subpopulations, weighted by the sizes of the subpopulations, the variance of the allelic state in the total population. Which is not provided by amova. Thus you cannot compute the fixation index from amova results.


How are UMIs used to dedupulicate in Drop-seq tools? The module DigitalExpression which is part of the popular Drop-seq tools digitally count gene transcripts.

The manual is not very clear on how exactly it resolves the UMI information.
Maybe I'm missing something but I was wondering if UMI are used in isolation or in conjunction with the barcode and or transcript?

The manual i.e. Drop-seq Alignment Cookbook only states:

To digitally count gene transcripts, a list of UMIs in each gene, within each cell, is assembled, and UMIs
within edit distance = 1 are merged together. The total number of unique UMI sequences is counted,
and this number is reported as the number of transcripts of that gene for a given cell.

To clarify my point, I'll use an example. Say my reads are only 4 bp and read 1 has the cell barcode at position 1+2 and UMI at position 3+4.  read 2 constitutes the transcript read.

Given the following example reads:

seqence1:   read1 = AATT  read2 = ATGC
seqence2:   read1 = AATT  read2 = TTTT

UMI and barcode are the same but read2 is different, would they be collapsed to 1 UMI count or 2?

Similarly as case with same UMI and transcript but different cell barcode:

seqence1:   read1 = AATT  read2 = ATGC
seqence2:   read1 = GGTT  read2 = ATGC

Also of note, UMI deduplication is the default setting as of the parameter below:

OUTPUT_READS_INSTEAD=Boolean  Output number of reads instead of number of unique molecular barcodes.  Default value:
                              false. This option can be set to 'null' to clear the default value. Possible values:
                              {true, false}

UMI and barcode are the same but read2 is different, would they be
  collapsed to 1 UMI count or 2?

If the two reads share cell barcodes and UMI sequences and align to the same gene, (even if they are different sequences) then they originate from the same RNA molecule, and should only count as 1.

https://support.10xgenomics.com/single-cell-gene-expression/library-prep/doc/user-guide-chromium-single-cell-3-reagent-kits-user-guide-v3-chemistry


What does "motif 2 bp or longer" mean? I'm reading a preprint on a new somatic variant discovery tool, and the first sentence of the Discussion section has me confused.

  Across the four datasets analyzed in this study, we discovered that the major source of disagreement between callers originates from somatic variants called within STRs, in particular if the motif is two base pairs or longer.

The "two base pairs or longer" bit is what confuses me. I interpreted "motif" as referring to the short tandem repeats (STRs), which are by definition 2 bp or longer. A motif of < 2 bp is...a homopolymer, not really a motif in the common sense.

Am I missing something here? Or is this simply an awkward word choice on the part of the authors?

What's your reference for the definition of STRs? I think it is still ambiguous among the community. Wikipedia states motifs of >=2 bp. However other references include homopolymers as well:
Source 1, source 2 and source 3.


clusterProfiler, groupGO: meaning of GeneRatio I am running groupGO function from clusterProfiler:

ggo <- groupGO(gene     = head(gene, 200),
           OrgDb    = org.Mm.eg.db,
           keyType = "ENSEMBL",
           ont      = "BP",
           level    = 2,
           readable = TRUE)

And getting the following output:

> head(ggo)
                   ID            Description Count GeneRatio
GO:0000003 GO:0000003           reproduction     5     5/200
GO:0008152 GO:0008152      metabolic process    40    40/200
GO:0001906 GO:0001906           cell killing     0     0/200
GO:0002376 GO:0002376  immune system process    20    20/200
GO:0006791 GO:0006791     sulfur utilization     0     0/200
GO:0006794 GO:0006794 phosphorus utilization     0     0/200

As can be seen cell killing, sulfur utilization and phosphorus utilization have zero genes. Then, why are they here at all? Should not the groupGO show only biological processes the fed in genes are related to? This seems very weird to me. What is the meaning of the GeneRatio column?

As you seem to be learning, the clusterProfiler package is a bit "special". It doesn't filter the results you're seeing for relevance, it returns everything. The GeneRatio is the number of genes in your input list associated with the given GO term (e.g., 5) and then a / and the total number of input genes. This is generally not that useful. Since many things aren't documented, you may find the source code helpful.


Viral genome finishing I have assembled poxvirus genome using Ray. The assembly is good. Out of several thousand contigs I got, I was able to get one scaffold using Contiguator tool, which is about 90% of my genome. I have well annotated reference sequence. If there any good tool for filling gaps and checking assembly quality?

Things I've tried so far:

IMAGE from PAGIT. It asks for ace file which I guess is obtained from Newbler. So, I am not able to use this.
gapfill_py script from broad: There are some gaps still missing.

I came to know that Consed is good. But, I have windows which is not compatible for Consed.

I checked my BAM file with Tablet. There are reads that cover almost entire gene. Are there any good tools or manual ways to finish my assembly so that I can get the full genome?

Assuming that gap filling is a domain-independent and works the same for viruses as for bacteria or eukarytes.

There are two real options I am aware of: Sealer a successor of GapFiller from ABySS package and GapCloser from SOAPdenovo2 package. According to Sealer paper it's better, but I am not sure if any indepented benchmarking have confirmed their conclusion. On my insect data GapCloser worked better than GapFiller (back then I have not tried to compare it to Sealer).


how to extract a different substring for each row of a dataframe in R? I have a dataframe with putative alignments of thousands of probe sequences. I would like to pull out the adjacent nucleotide from the genome for each alignment. My dataframe includes seq.name for the target chromosome and I've added an adjacent_pos variable with the location on the target chromosome that I want. Using read.fasta from the phylotools package gives me a genome in the format of a dataframe with seq.name and seq.text for each chromosome and unplaced scaffold.
Working on a test data set with a single chromosome, I successfully used something like this:
adj <- cbind(alignments$adjacent_pos, alignments$adjacent_pos)
alignments$adjacent_base <- str_sub(genome$seq.text, adj)

Now back to the whole genome, I can't wrap my head around how to direct str_sub to the correct seq.name for each seq.position. I'm sure there's an elegant solution, but I haven't been able to find it!

Say your file looks like this:
cat test.fa
>chr01
ACCGGTT
>chr02
ATGCATGC

Something like what you have done, to read in the fasta and extract the sequence. Key part is creating a GenomicRanges object:
library(GenomicRanges)
library(BSgenome)
library(Biostrings)

adj <- data.frame(chr=c("chr01","chr02"),start=c(2,3),end=c(2,3))

    chr start end
1 chr01     2   2
2 chr02     3   3

adj = makeGRangesFromDataFrame(adj)

fasta <- readDNAStringSet("test.fa")
SEQ = getSeq(fasta,makeGRangesFromDataFrame(adj))

as.character(SEQ)
[1] "C" "G"

or if the fasta file is huge, you can use Rsamtools, indexing the fasta file first:
library(Rsamtools)
fasta = FaFile("test.fa")
indexFa('test.fa')

as.character(getSeq(FaFile("test.fa"),adj))
chr01 chr02
  "C"   "G"