Getting the licenses of all conda installed packages

To my knowledge there's no single listing of the licenses, but you can:

conda list | awk '{if(NR>3) printf("%s=%s", $1, $2)}' | xargs conda info | grep license

Allegedly conda info --license should work, but at least in conda 4.6.3 it doesn't show any license information. It unfortunate that there's no convenient way to get a tabular list of package:version:channel:license information from the command line.


Extract reads from bam files by their @RG

Typically I use samtools for operations like this.  Specifically I use samtools view with either -r or -R flag depending on the use case.

-r STR
Output alignments in read group STR [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.

-R FILE
Output alignments in read groups listed in FILE [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.


Finding the nearest gene to a specified gene region

You're looking for bedtools closest, which thankfully defaults both considering overlapping genes and outputting all hits in case of ties.

  bedtools closest -a foo.txt -b genes.bed | awk '{foo = sprintf("%s:%i-%i", $6, $7, $8); OFS="\t"; print $1, $2, $3, $4, $5, foo}'

You'll need to strip the header and ensure the file is sorted. I've piped the output to awk to reformat things. The output then is something like:

chr1    24613189    24613641    0   3   chr1:24613188-24613971
chr1    24613971    24614812    0   4   chr1:24613973-24614651
chr1    35788729    35789387    -5  0   chr1:35806973-35810462
chr1    93151586    93160225    0   -2  chr1:93151348-93160948
chr1    93151586    93160225    0   -2  chr1:93151582-93160870
chr1    107511454   107525597   -2  0   chr1:107511422-107525598

Those results are for the Gencode M15 release (the BED file I used has transcripts rather than genes, so you'll get a bit different coordinates). As an aside, the gene name would normally be output as well, just add column 9 to the output.


Is there a software or web server to calculate thermodynamic parameter enthalpy change of cellulase?

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


How can I find mutations associated with disease in human histone residues?

Ensembl contains this information: When you go to the “phenotype” menu item of a given gene, you will see a list of variants (potentially after clicking on “ALL associated variants”) with their associated genomic position.

You can subtract the transcript’s start position from that position to find out the residue (of the translated sequence, and be careful to subtract the lengths of all preceding introns!). Unfortunately this seems to require some manual work.

Alternatively you can query the individual SNPs, which will tell you their position in all transcripts (in this case: Leu 298/68/244):


In SBOL, when should I annotate a DNA sequence vs. making sub-components?

Typically you want to use a ComponentDefinition when an sub-structure is something you might want to pull out and re-use in another genetic design, and a SequenceAnnotation whenever it is not.

For example:

Promoters, terminators, and protein coding sequences are typically best represented with their own ComponentDefinition, since these are often pulled out and re-combined in new genetic designs.  Representing such an element with a ComponentDefinition makes it easy to find the commonalities between designs.
Assembly scars, binding sites, and cut sites are typically best represented with a SequenceAnnotation because they are typically not very interesting except in the context of their surrounding genetic context. Representing such an element with a SequenceAnnotation helps keep them associated tightly with that context.

Not every case is cut and dried, however.  For example, ribosome entry sites are often treated as separable components, but their performance is very tightly tied to the specifics of the coding sequence that they modulate.  Thus, in some cases it may make more sense to represent them separately as a ComponentDefinition and in others to represent them as a SequenceAnnotation on a joint RES/CDS ComponentDefinition.


python does not quit when the input file is too big

You've hit the RAM bottleneck. Its very common in de novo assemblies and resolved using template assemblies. I don't know this program but suspect it is performing a pairwise analysis for every piece of data in a genomics data set.
The solution is to either:

find a bigger machine, e.g. a big multi-core machine and make sure the program distributes across every core.

parallelise the calculation (an excellent way to overcome RAM problems), or,

find an approximation to the calculation that is less computationally demanding, templates are a good way to circumvent this.


Nvidia Parabrick fq2bam pipeline error - No such file or directory

My query was answered by NVIDIA Parabrick developer team. The installation instruction are available at https://forums.developer.nvidia.com/uploads/short-url/t5uXsZgdMXx0G2wDcYamFl9sO2x.pdf.


Filtering out reads from a reference (e.g. rRNA) using bowtie2

If you want to remove rRNA, use tools that are made for this. We usually go for SortMeRNA, which can do exactly what you describe your goal is, i.e. go from FASTQ to FASTQ, splitting your reads into rRNA and non-rRNA.
In our experience, it is very accurate, if a bit slow at times.


How to identify genes from a genome assembly of C. Elegans?

You don't want BLAST here, or at least not regular BLASTn. A better method would be to use tBLASTn to map the known C. Elegans proteins to the translated genome. However, I would recommend using a slightly more sophisticated approach:

Collect fasta sequences of all known C. Elegans proteins (amino acid sequences, not nucleotide). This should be easy enough to do on WormBase. For example, the protein sequences for C. elegans strain VC2010, assembly PRJEB28388, can be found here: ftp://ftp.wormbase.org/pub/wormbase/releases/current-production-release/species/c_elegans/PRJEB28388/c_elegans.PRJEB28388.WS279.protein.fa.gz

Use a tool that can model introns and exons to map a protein into a genome. I haven't worked in this field in more than 10 years now, so it is very likely that there are other programs around today, but back in the day, I would do this using either exonerate or genewise. With exonerate, the command would be:
exonerate -m protein2genome -n 1 -t assembly.fasta -q allproteins.pep > out

That will give you a list of best-hits for each of your input proteins and should accurately map the protein to the genome. Once you have that, you can extract the matches and start comparing.


How can I colour boxes in Gviz AnnotationTrack in R?

Did you try the fill argument?

Something like this:

track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'), fill=c('red','blue', 'red'))


RNA-Seq: clustering/treatment of genes with low expression

rlog(normalized counts) is going to be more robust than log2(TPM), so use it instead.
Do it afterward, keeping the low counts can be helpful for library-size normalization.
This is very much more of an art than a science, though have a look at the kOverA() function from the genefilter package. That will give you a bit more fine-grained control on filtering.
You'd be better off using something like zFPKMs if you want to filter by "expressed or not", since there's no really great definition of that.


How can I count the number of reads that support a variant in a bam file?

If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically:

samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam |
    cut -f 5 | tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c

That'll give you a count of how many of each base were seen.


RNA-seq: How to get new expression count after normalization

If I'm interested to get normalized gene counts, can I go ahead and
  multiple each individual's norm.factor into its gene counts? For
  example, the expected count for IND1 for Gene1 is 100 and it's
  norm.factor is 0.80, can I say that the normalized gene count is
  100*0.80=80?

No; you divide, but if you poke around, you can probably find a way to get edgeR to return the normalized counts.

  I think since normalization uses per person library size, each
  person's normalization is independent of the others.

Not true for the default normalization used in edgeR, TMM.  CPM normalization is independently determined for each sample.

In general, normalize everything together, unless you think that the algorithm assumptions might be violated by including all of them together.


Can I use my computer while MinION is sequencing without negatively affecting the run?

Yes, that should be fine, especially if basecalling is disabled. Locking a screen should be fine as well.

I prefer to do basecalling after a run anyway, because there's a chance that it could chew up too much processor time and result in lots of skipped reads.


Without control data from your subjects, I don't think there's really no way to distinguish somatic mutations from germ-line mutations. The best you can do is to screen out common variants, which are germ-line mutations that are shared by large numbers of individuals using the population frequencies from something like the  Exome Aggregation Consortium:

annotate_variation.pl -downdb -webfrom annovar -build hg19 exac03 humandb/
annotate_variation.pl -filter -build hg19 -dbtype exac03 example/ex1.avinput humandb/


What is a poly(T) primer anchor sequence?

As shown in figure 1 of that paper, the anchor sequence is additional sequence that is added to the end of the polyT primer that provides a known and non-target binding site for additional PCR reactions. I gave a presentation about a month ago about the process that Oxford Nanopore Technology use in their strand-switching cDNA kit, which may be useful for understanding the process used by Tang et al.. Here are my strand-switch/PCR images from that presentation:

PolyT primer addition, with mismatched anchor sequence

Strand switch primer addition, with overhanging anchor sequence

PCR amplification, using primers that are complementary to the anchor sequences

In the case of the Tang et al. paper, they use a polyA-tailing method for strand-switching, but the process is essentially the same.


The start position is the first correctly mapped base. One way to do it, is to use GenomicAlignments in R to read in the alignments.

For example the sam file looks like this:

cat test.sam
@HD VN:1.3  SO:coordinate
@SQ SN:chrI LN:228861
D00689:415:CD4LUANXX:3:1316:7447:98083#CGGCTATGCAGGACGT 0   chrI    769 0   42M8S   *   0   0   CCATAGCCCGGTGGGGCCTGTGCTGCCACAGAGGGTGTGTCTAGATCGGA  B/B//FFFFFFBFFFF//7FF/<F/FFBFFFF//<F///<F/BFB//<F/  NM:i:1  MD:Z:21A20  AS:i:37 XS:i:37
D00689:415:CD4LUANXX:3:2308:9597:30445#CGGCTATGCAGGACGT 0   chrI    3965    0   3S47M   *   0   0   CCACACACCCACACACCACACCCACACCACACCCACACCACACCCACACC  BBBBBFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFFFFBF<FFFFF<FF  NM:i:0  MD:Z:47 AS:i:47 XS:i:47

We convert to bam:

samtools view -bS test.sam > test.bam

And read using R:

library(GenomicAlignments)
gr = readGAlignments("test.bam")
gr

GAlignments object with 2 alignments and 0 metadata columns:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     chrI      +       42M8S        50       769       810        42
  [2]     chrI      +       3S47M        50      3965      4011        47
          njunc
      <integer>
  [1]         0
  [2]         0
  -------
  seqinfo: 1 sequence from an unspecified genome

The start and the end can be retrieved using:

start(gr)
[1]  769 3965
end(gr)
[1]  810 4011


--thread 4 is essential, --maxiterate just try and get it working, not getting it to work optimally. Just make sure there are no other processes running on the machine when --thread 4 is operating, e.g. you are using another machine as your desktop
The gist is as you spread the calculation across your processors (4) and simplify the calculation (lowering the iterations), you are getting further into the calculation, for example now you are getting into gap penalties.
If you are agreed then I would try performing the alignment using just one chromosome of the human genome against the chimp genome and I suspect that would work. Then take it from there.
The objective is to get the calculation working, for example I wouldn't trust a result with max. iterations of 2, but once you have it working you can then start refining the parameters towards optimal results. The i3 processor isn't perfect in context.
Using --thread 4 is super important and you want to increase maxiterate towards the default if the machine supports it.
I suspect the solution is to perform the calculation on an 8+-core machine and set --thread 8 (or higher depending on the number of cores on the machine), at least it will get further into the calculation. With 4-cores and an i3 (diagnostics of RAM usage would be nice) in any case it just doesn't seem to like it. The alternative is to perform this across the genome via sampling the chromosomes on your current i3.

Just to explain why this works. When you have a RAM intensive job and you spread it across cores the RAM overhead massively reduces for each core that is added. Don't ask me how this works in the hardware mecahanics of RAM distribution per core. The reason I know this is because of emperically parallelising Bayesian calculation using each chain on a different core.


Why don't the "child" peptide masses sum to the parent in this De novo sequencing example?

In a peptide tandem mass spectrometry experiment, each fragment (child) ion is typically the result of breaking a single bond in the peptide backbone yield b and y ions. Because each ion fragment ion is the result of one fragmentation event, we don't expect the fragment ion masses to sum to the precursor mass.
If your peptide was fragmented well, you should get a ladder of b and y ions where the differences in masses between ions in the b and y ion series correspond to the masses of amino acids (+/- any modifications and neutral losses, such as the NH3 and H2O mentioned above).
A great resource for understanding these concepts is this tutorial:
https://www.ionsource.com/tutorial/DeNovo/introduction.htm


how to extract a different substring for each row of a dataframe in R?

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"


RNASeq: Normalization, stabilization, gene length and rlog

I am not sure what you mean by "decimal data", however, if you really want to use "gene length" information for normalization, take a look at EDASeq. Unlike DESeq or EdgeR, its normalization step takes gene length into account to produce "normalized counts", which you can feed into DESeq or EdgeR according to their manual.

Below is from the said manual:

  The normalized counts (or the original counts and the offset) obtained
  using the EDASeq package can be supplied to packages such as edgeR
  (Robinson, McCarthy, and Smyth 2010) or DESeq (Anders and Huber 2010)
  to find differentially expressed genes.


Detecting SARS-CoV-2 variants from the mixed virus population

Based on a quick request from someone else, I wrote some code a couple of years ago for assigning read groups based on variants at particular locations; this was set up for HIV variant splitting, so I expect it should also work on SARS-CoV-2 nanopore reads as well:
https://gitlab.com/gringer/bioinfscripts/-/blob/master/samVarSplitter.pl
It processes a SAM stream [i.e. in text mode] using Perl, so might be a bit slow, but should be fine for a few hundred thousand sequences. You'd use it as follows:
samtools view -h mapped_reads.bam | \
  samVarSplitter.pl [-ref <refName>] [-pos <int>] | \
  samtools sort > grouped_reads.bam

Then the reads can be split into different files based on their group using samtools split:
Usage: samtools split [-u <unaccounted.bam>] [-h <unaccounted_header.sam>]
                      [-f <format_string>] [-v] <merged.bam>


Searching tool to calculate phase/switch error rate

Assuming your files are as you show, and the genotype field is always the 1st entry of the sample information field (the 10th), then you can just do this directly with a simple awk script:

$ awk -F"\t" -vOFS="\t" 'BEGIN{print "#CHROM","POS","ID","REF","ALT","Truth","Test"}{ if(/^[^#]/){var=$1$2$4$5; gt=substr($10,1,3); if(NR==FNR){a[var]=gt; next;} if(var in a && a[var]!=gt){print $1,$2,$3,$4,$5,gt,a[var]}}}' truth.vcf test.vcf
#CHROM  POS     ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0

Here's the script split into separate lines for clarity (you can still copy paste directly into a terminal):

awk -F"\t" -vOFS="\t" '
    BEGIN{
        print "#CHROM","POS","ID","REF","ALT","Truth","Test"
    }
    {
        if(/^[^#]/){
            var=$1$2$4$5;
            gt=substr($10,1,3);
            if(NR==FNR){
                a[var]=gt;
                next;
            }
            if(var in a && a[var]!=gt){
                print $1,$2,$3,$4,$5,gt,a[var]
            }
        }
    }' truth.vcf test.vcf

The basic idea is:

awk -F"\t" -vOFS="\t": set the input and output field separators to a tab.
BEGIN{print "#CHROM ... }: print a header.
if(/^[^#]/){: skip the header lines of the input vcf files
var=$1$2$4$5: get a "signature" for this variant, by concatenating the chromosome, position, ref and alt fields.
gt=substr($10,1,3); : get the genotype value.
if(NR==FNR){: if this is the 1st file/
a[var]=gt; next: save the genotype for this variant and skip to the next line.
if(var in a && a[var]!=gt){: this part will only be run for the 2nd file because of the 'next' in the previous block. If this variant was in the truth set and if its genotype there isn't the same as in the test set, print the variant and the two genotypes.

If you need the error rate, you can extend the script thusly:

awk -F"\t" -vOFS="\t" '
BEGIN{
    print "#CHROM","POS","ID","REF","ALT","Truth","Test"
}
{
    if(/^[^#]/){
        var=$1$2$4$5;
        gt=substr($10,1,3);
        if(NR==FNR){
            a[var]=gt;
            next;
        }
        tot++;
        if(var in a && a[var]!=gt){
            bad++;
            print $1,$2,$3,$4,$5,gt,a[var]
        }
    }
}
END{
    printf "Correct: %s of %s (%s%)\n",tot-bad,tot,(tot-bad)*100/tot;
    printf "Wrong: %s of %s (%s%)\n",bad,tot,bad*100/tot;
}' truth.vcf test.vcf

On your example dataset, that prints:

#CHROM  POS ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0
Correct: 0 of 5 (0%)
Wrong: 5 of 5 (100%)


How to combine multiple files into one file?

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 obtain only one cluster's data (.csv file) in Seurat

You can access data within the Seurat object using GetAssayData, and extract a list of cell names for the cluster you're interested in using WhichCells:

library(Seurat)

cells.use <- WhichCells(object = pbmc_small, ident = 1)
expr <- GetAssayData(object = pbmc_small, assay.type = "RNA", slot = "data")[, cells.use]
expr <- as(Class = 'matrix', object = expr)
write.csv(x = expr, file = "expression_cluster1.csv", quote = FALSE)


I found out that the package PyPDB has the function get_ligands that does exactly what I'm looking for.


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

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


How can I add for several bed files the header : track type=narrowPeak name=“narrowPeak” preferably in python ,can handle with R

I don't think you need to go to R or Python to do that. In bash
set -euxo pipefail
for f in *.bed
do
  bn="$(basename ${f})"
  on="${bn%%.bed}_with_header.bed"
  echo 'track type=narrowPeak name="narrowPeak"' > "${on}"
  cat "${f}" >> "${on}"
done

EDIT: I just read you want to set the name based on the filename.
EDIT2: And also to add an arbitrary directory
set -euxo pipefail

INDIR="/path/to/bedfiles"

for f in $(find "${INDIR}" -name "*.bed")
do
  bn=$(basename "${f}")
  on="${bn%%.bed}_with_header.bed"
  dn=$(dirname "${f}")
  echo "track type=narrowPeak name=\"${bn}\"" > "${dn}"/"${on}"
  cat "${f}" >> "${dn}"/"${on}"
done


Compare multiple alignment results' aligned bases

I don't know of a single tool that can read multiple sets of BAM files and compare them. But you can use the Sam_Reader python3 module to create a csv file that summarizes each position in the genome for each of your 3 results and then compare them manually.

It should be installable with:

pip3 install sam_reader

And then generate the CSV with

from sam_reader import Sam_Reader
my_files = Sam_Reader('my_bam_file_or_folder/')
my_files.per_base_stats(write_file='new_file.csv')

*but FYI this module requires pysam, which does not work on windows


Both of the annotation tools you used are designed for prokaryotic genome annotation. I would not expect them to work very well for viruses.


Retain only part of a file name and fasta header in fasta directory

I think you're working too hard to capture the uce: expression.
Try
sed -Ee '1s/^.+uce:(uce-[0-9]+).+/\1.unaligned.fasta/' file_glob

to replace the first line content. You may want to add the in-place flag after testing the expression on a few files.
For the mv, use command substitution, which runs the command in a sub-shell and uses the output value as the parameter:
for f in file_glob
do
    mv -- "$f" "$(echo \"$f\" | sed '...')"
done


Correlating gene expression with qualitative variables

Logistic regression would generally be a bad choice for non-binary outcomes. In such cases, linear regression (or a GLM more generally) still works fine. You can already do that in the standard R packages for RNAseq (DESeq2, edgeR, and limma), where the fold-change is in whatever units you're measuring your quantitative trait in.


Determine if an organism is microbic or unicellular

I would start by saying that your statement that no animals/plants are microbes is false. There are indeed microbial algae (plants) and microbial animals.

So it is indeed difficult. (I'll briefly note that unless your protocol specifically gets rid of environmental DNA, you will probably have some contaminating non-microbial DNA no matter what.)

However, we can use a different version of your proposition to make a simple (over-simple, most would argue) approximation: all prokaryotic organisms, e.g. bacteria and archaea, are microbial, and all eukaryotic organisms are not. This is not hard to do if you are already using NCBI taxonomy, as you can use e.g. this mapping. It has the advantage that all bacteria and archaea are definitely microbial.

Of course it's quite wrong, as you already know. In principle you could set up a complex decision tree sorting this out, though I imagine that's the thing you wish to avoid.

In terms of other genome databases, you might consider looking into GOLD, which has extensive metadata in a slightly less arcane setup than NCBI, which might be helpful. But there is a high rate of missingness there, and also there is still not necessarily a great "microbe"/"not microbe" mapping.

There are also existing databases of microorganisms, which may or may not satisfy your requirements; I have little familiarity with most of those so caveat emptor.

I believe (hope) that NCBI taxid is a good way to do lookups in all such databases, but it may not be true.


Denovo, Stacks: Getting an “ambiguous redirect” error

This most likely means your variable $logfile does not exist.


Methods to predict one connected structure from two known separate structures

There are two approaches here.
Physiological
If domain A and B of your protein interact with each other physiologically, then you need to protein-protein docking. This is not too straightforward as protein are flexible.
Several online servers exists such as Haddock and include ClusPro, which use more traditional approaches. While the online server DEMO (I-Tasser group) does a structural homology based search and a final MD based repacking to better minimise the structures (no protein-protein docking).
Many of these work by magic as far as the user is concerned. But understanding what is going on is really important, even if it is overly laborious. When you don't know where a protein goes, it's called blind or global docking, which is computational expensive and so imprecise that the effort required to prove it is correct for peer-review is too steep for it to be useful. If you are doing local docking it is okay, which may be the case here depending on the linker. Specifically, you need to generate more 'decoys' (potential poses) the more unsure you are: tens of thousands or more. Then you have find clusters and see how each scores following more precise refinement. As an example of a tediously and manually set up run, Rosetta has several docking protocols, which explain the steps required —of note are the many judgement calls, such as in the local refinement, where one has to judge whether setting -dock_pert 3 8 (3 Å and 8° perturbation) is too strict etc.
Pretty picture
The protein-protein docking approach will result in a globular protein, which may be incorrect and majorly is not great for pictures.
If you want a pretty picture of two or more domains most likely you would want them placed like beads on a string. If you have a linker that is not solved the recommended distance to space them is $(3.5 N_{missing})^{0.5}$ Å —unfortunately I do not know where this lore comes from. If it helps I wrote a script that does this in automatic but manual placement works best.
Linkers
In the situation you have a good confidence about how the two parts fit together (or have placed them like beads on a string), you may want to add linkers. If you absolutely don't care about correct bond topology here is a blog post on how to join bits (shoddily) with PyMOL and here is one to do it properly with Rosetta or Pyrosetta.
If one wanted to do a better job, one tool to use is Rosetta remodel, which has a bit of a learning curve, but is very powerful for protein design. This requires the PDB to be on a single chain (pymol alter command, followed by create will do the trick: more here) and actually the amino acids to not have gaps —PDB number being the same as the pose numbering. Rosetta gives better results if the poses are energy minimised, but for this it does not matter too much. Then a "blueprint file" formatted as discussed in the documentation. Briefly, in a way that residues left alone are in the form 1 M . (where residue 1 is a methionine), while the residues before and after the insertion are like 100 D L PIKAA D (meaning D100 is changed to a Asp in loop SS; not actual numbers or residue), while between entries are added 0 X L PIKAA G etc. To model loops, -generic_aa G is a good option as this tells the default rough fitting (centroid mode) to use glycine and not valine as generic residue.
The top pose will actually not have a superb score, so a 3 (regular) or more cycles (15 thorough) with Rosetta Relax application either of the whole protein or with a movemap specific to the linker, may be done.
There actually is an option in Rosetta remodel to have a domain insertion, i.e. it would swing around like a medieval flail. However, it will not repack and therefore properly dock so this option is a bad idea for this specific application.


How to download gene expression data from NCBI gene database

You can download raw data from SRA (NCBI). On this website under RNA-Seq alignments, you'll find the samples. Click on them, for example SAMN05231885, and then click on PRJNA325427 to see the link to the SRA data (click on 1). Here you can download the raw fastq files of the experiment.


Error in adding data in columns of R data frame

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

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


How to import large .bed, .gff, .vcf, .paf, .sam files into an SQL database?

Answer from @liam-mcintyre converted from comment:
I don't use dask as it doesn't support enough pandas functionality (unfortunately). With pandas I do it with read_csv... if its big then read in chunks and send chunks to separate threads. If you want to ask a specific question with example data etc then I can show code.


How to subset a VCF by chromosome and keep the header?

This seems to work for me with bcftools filter and the -r or -R argument.

   -r, --regions chr|chr:pos|chr:from-to|chr:from-[,...]

  Comma-separated list of regions, see also -R, --regions-file. Note that -r cannot be used in combination with -R.

   -R, --regions-file FILE

  Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The columns of the
             tab-delimited file are:
             CHROM, POS, and, optionally, POS_TO, where positions are 1-based and inclusive. Uncompressed files are stored in memory, while
  bgzip-compressed and
             tabix-indexed region files are streamed. Note that sequence names must match exactly, "chr20" is not the same as "20". Also note
  that chromosome
             ordering in FILE will be respected, the VCF will be processed in the order in which chromosomes first appear in FILE.
  However, within chromosomes, the
             VCF will always be processed in ascending genomic coordinate order no matter what order they appear in FILE. Note that
  overlapping regions in FILE can
             result in duplicated out of order positions in the output. This option requires indexed VCF/BCF files. Note that -R cannot be
  used in combination with -r.

Here's an example:

$ bcftools filter vcf_nocomp_merge_geno98.vcf.gz -r 4 | head -n 38 | colrm 100 1000000
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20171222
##source=PLINKv1.90
##contig=<ID=1,length=249167691>
##contig=<ID=2,length=242695901>
##contig=<ID=3,length=197800245>
##contig=<ID=4,length=190915651>
##contig=<ID=5,length=180666277>
##contig=<ID=6,length=170877445>
##contig=<ID=7,length=159086805>
##contig=<ID=8,length=146293415>
##contig=<ID=9,length=141018424>
##contig=<ID=10,length=135434552>
##contig=<ID=11,length=134938471>
##contig=<ID=12,length=133763353>
##contig=<ID=13,length=115045730>
##contig=<ID=14,length=107285438>
##contig=<ID=15,length=102369712>
##contig=<ID=16,length=90141356>
##contig=<ID=17,length=81006630>
##contig=<ID=18,length=78014583>
##contig=<ID=19,length=59071322>
##contig=<ID=20,length=62906515>
##contig=<ID=21,length=48077813>
##contig=<ID=22,length=51156934>
##contig=<ID=23,length=154847490>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on rea
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_filterVersion=1.2-187-g1a55e45+htslib-1.2.1-256-ga356746
##bcftools_filterCommand=filter -r 4 vcf_nocomp_merge_geno98.vcf.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  1_125   3_5555  4_7565  5_9
4   71566   rs13125929  T   C   .   PASS    PR  GT  0/1 0/1 0/0
4   87113   rs2006748   T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   110646  rs11727494  C   T   .   PASS    PR  GT  0/1 0/0 1/1
4   142550  rs11735742  T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   200631  rs6826124   A   G   .   PASS    PR  GT  0/1 1/1 0/1
4   221623  rs7695945   G   T   .   PASS    PR  GT  0/0 0/0 0/0


How can I find out the code DRIMSeq uses to make dmplotProportions objects?

I was curious how they have done it as well, so here expanded @Liopis comment :

wget http://bioconductor.org/packages/release/bioc/src/contrib/DRIMSeq_1.6.0.tar.gz
tar -zxf DRIMSeq_1.6.0.tar.gz
# if you need to find it, or use grep
# find DRIMSeq -name "dm_plotProportions.R"
less DRIMSeq/R/dm_plotProportions.R


R is a reference to R0, otherwise known as the basic reproductive rate, and means the number of new cases from a single patient infection. R0 has to be above 1 for a disease to persist. It is a classic equation which account for the probability of transmission against immunity. The y-axis is the number of infecteds, very poorly presented.

The idea is that R0 will be lower if there is containment (less people to infect because you can't move around), but will last longer, possibly because the number of naive patients (not encountered the virus) will be larger over a longer period of time.

That particular graph appears to be saying that everyone at risk of infection will get infected with or without control measures (the surface area under the graphs is approximately the same*) and the only thing that will halt the virus is acquired immunity. If this is mainstream economic thinking then I beg to differ, albeit I note the source is CDC (again see note*).

For example, all other respiratory viral infections are seasonal. If COVID-19  is seasonal containment is the right choice, because the window of transmission is limited.

*, Given this is of CDC origin I think either its a mistake, or at the very least was taken out of context, because the idea is that overall there is a lower number of people infected as a result of "social distancing" and that does not appear to be what the two curves show.


RNA_Seq data aligned used uniquely or multi mapped reads impact on result interpretation

There are several methods/tools that take into account multi-mappers and deal with them in reasonable ways (e.g. pseudo-alignment methods like salmon and more traditional alignment-based methods like RSEM). If you completely ignore multi-mappers, then yes, you will be losing information, which may or may not be valuable to you.
Both of those methods utilize the expectation-maximization algorithm (or a variant thereof) to determine the relative contributions of multi-mapping reads to the various mapping locations. With salmon, you can forego alignment completely and end up with very accurate gene quantifications in hand in like 20 minutes if you want.


GATK CombineVariants complains the contig order in the VCF files

Answered in the GATK forum

  Try regenerating the index files for your VCFs. Picard SortVcf doesn't do it for you iirc, so when GATK looks at the index files (before opening the VCFs themselves) it sees the old order and complains.


How can I take cell number into account to find total RNA expression?

(I currently cant comment as I used up all my rep for a bounty.)
Could you explain how you calculate the average? In you plot the averages are between an expression level of 2 and 3, while you report values over 10. Given the data you show, there is no significant difference and I would guess p-values are above 0.8!
And you mention that the violin plots suggest the control groups have higher expression.. but how do you see that? For me they are identical with minimal trends that the PLX5622 groups are higher.
And to address your question regarding showing the total: No, you would want to show the mean or median, as the total is pointless, given deviating sample sizes.
EDIT: Overall, I have a feeling you uploaded the wrong picture!
EDIT2: I think I got your point now: The difference in mean is caused by an increased number of outliers in your controls at the value 0. Despite those, the means are identical and this is what can be seen through the violin plots. You should rather try to diminish the effect of these outliers.


Comparing groups of 16S rRNA sequences - how to?

If you want to make an all vs all with blast this can be helpful. You can create a blast database with your sequences and then search against itself.
makeblastdb -in 16S_sequences.fasta -dbtype nucl -out my_16S_sequences_db

blastn -db my_16S_sequences_db -query 16S_sequences.fasta -outfmt 6 \
       -out allvsall_16S.tsv -num_threads "$(nproc)"

If you want a percent identity matrix for a pairwise alignment you can try with this approach: https://www.biostars.org/p/220154/


The coefficient is the third column from the right, you can tell because some are positive and other negative. When the value is negative the association is a negative correlation. When it is positive (last two) that is more like the result you are seeking. The further away from zero the stronger the correlation (positive or negative)
These are very large values for the coefficients. The SNP is the last column on the right. I would disgard the top two results because the probabilities are too borderline, e.g. if you used Bonferoni's correction.
There isn't enough information on this calculation nor the data set to make any further insights other than you have 7 SNPs showing PRS


Staff from rcsb kindly reminds me the solution - add the following lines in the json query.
"request_options": {
"return_all_hits": true
}


How do you write a .gz fastq file with Biopython?

I'm not sure I'm doing it the best way, but here is an example where I read a compressed gzip fastq file and write the records in block gzip fastq:

from Bio import SeqIO, bgzf
# Used to convert the fastq stream into a file handle
from io import StringIO
from gzip import open as gzopen

records = SeqIO.parse(
    # There is actually simpler (thanks @peterjc)
    # StringIO(gzopen("random_10.fastq.gz").read().decode("utf-8")),
    gzopen("random_10.fastq.gz", "rt"),
    format="fastq")

with bgzf.BgzfWriter("test.fastq.bgz", "wb") as outgz:
    SeqIO.write(sequences=records, handle=outgz, format="fastq")


In a tissue or cell line that used to be diploid and has undergone genome wide doubling, you expect more variants having observed allele frequencies that would also occur in tetraploid organism. That is, observed VAF distribution will have (possibly small) peaks at 0.25 and 0.75, in addition to the diploid peaks at 0, 0.5 and 1. Compare the distributions in the tissue or cell line of interest to those of the control, that is, known diploid. Also if it is a doubled male genome, the sex chromosomes will show a new peak at AF 0.5, that is, heterozygous X or Y variants. This is just an idea, I have never tried to look for such events.


You can try Mauve multiple genome aligner.


What are Approximate Read Counts (Library Sizes) and Lengths (Insert Sizes) for Next-Generation DNA Sequencers?

Both of those depend on the application of the library.  For example whole genome assembly of a mammal would require a larger library size than targeted re-sequencing of a virus.  Also insert size varies if you are doing genome assembly a larger insert size is valuable to allow for long range linking and scaffolding.

I think that if you edit your question of be more specific people will be able to provided a more concrete answer other wise it is "it depends"


How to adapt the fgseaL function to perform rapidGSEA computation of gene ranks across 9 different phenotype labels?

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

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

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

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


Error in PSI-BLAST when running psipred locally

This error could be due to the fact that you are using legacy blast, have you tried using BLAST+ instead?
If you wish to use PSIPRED with BLAST+, then use the runpsipredplus script
here, rather than the normal runpsipred script.


How to print alternative allele from VCF with htslib c++

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


Is there a single cell RNAseq equivalent of GTEx or TCGA?

I think the most well-known scRNA dataset that is similar to GTEx in the Human Cell Atlas (https://www.humancellatlas.org). It is pretty much still ongoing but there are some preliminary data.


Compare two networks of the same genes between two species

So, you can give a try to this Cytoscape app that will allow to perform comparison between networks.

This app is not specifically designed for comparing networks across different species but it will enable you comparing any kind of networks.

These are the papers they published, which can give you some hints.

About your specific problem...I can't help you about what to look for when you're comparing networks from different species (not an expert in this specific topic) but as far I know you should look for something related to patterns occurring across different networks (i.e. species), network motifs, distribution of network centralities, conserved structures. In general, to compare networks, I would look for global properties (e.g. diameter, clustering coefficient...) while I would use local properties (e.g. degree, closeness, and network centralities in general) to compare the role of a node that you retrieved in all your networks.

Also, since you're looking for some literature, these links may help:

About conserved patterns of protein interaction in multiple species
About comparing genomic expression patterns across species
This one is slightly out of topic but can give you
some other ideas about comparing networks.


Database of RNA-Seq and SNP genotyping data

Human RNA-seq data, like genetic data is very sensitive data as it can be used to identify people (and their relatives), so people will likely not distribute it without entering a formal agreement that ensures that people can be held accountable if the data is misused, stolen or made public. Depending on the consent form signed by the patients (or their relatives), they may not be allowed to share it at all. One possibility would be to start a collaboration with a group that has access to such data and ask them to run your analyses on their data. You can then have a shared publication and it's a win for everybody.


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

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

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

df1.it = NCList(df1.ir)

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

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

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

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

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

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


R's nnet package only supports fully connected neural networks with one hidden layer. This is the most primitive type of network. I doubt it will work well for promotors finding. In addition, such network won't give you useful interpretations.

If you want to explore neural networks, you should use a one-dimension convolution layer as is described in this paper. This layer effectively represents a positional weight matrix (PWM). You can know which motifs are used. With backtracking, you can also identify the coordinates of motifs. To deploy such a model, you need to be fairly familiar with deep learning and to learn one deep learning framework. For framework, you may start with keras. Once you get used to keras, you can implement your full promotor finder in less than 100 lines of python code. Alternatively, you may try dragonn. It is supposed to simplify deploying models for DNA sequences. I have no experience with it, though.

For promotor finding, it is also worth trying traditional methods. It is interesting that few/no neuralNet-based publications have evaluated traditional methods, probably because many ML people know little about classical motif finding. These methods could work well, too.


How do I find the most similar sequences from a large set of short sequences?

130k * 20bp is a small data set. At this scale, SSE2 Smith-Waterman may work well:

git clone https://github.com/attractivechaos/klib
cd klib && gcc -O2 -D_KSW_MAIN ksw.c -o ksw -lz
./ksw -a1 -b1 -q1 -r1 -t14 20bp.fa 20bp.fa > out.tsv

I simulated 130k 20bp reads from E. coli K-12. The last command line takes about ~2 CPU hours based on the partial output.

It is possible to make this much faster by using 5bp seeds. You can also try ssearch and fasta from the fasta aligner package. They are slower, probably because they are not optimized for such input. Computing Hamming distance with SSE2 will be much faster, but is probably of little use as it does not allow gaps. There are also Gene Myers' edit-distance based algorithms, which can be faster than SSE2-SW.

Bowtie2 (as well as most short-read mappers) won't work well at the default setting. Bowtie2 uses long exact seeds. It will miss many 1-mismatch hits, let alone 3-mismatch ones. The number you get from bowtie2 is an underestimate. You might be able to tune bowtie2, but to find 3-mismatch hits, you have to make it a lot slower.


Relationship between the instantaneous rate of change of a continuous time Markov model and the rate of a global clock model?

Essentially you are comparing a "local clock" with a global clock and whether the instance rate relates to the global it depends on the model and in MCMC clocks in most cases the answer is no. A global clock is a very specific and very limiting model (below). Everything below is a comparison between taxa not between alignment positions.

"Local clock" in the context of MCMC clocks is in inverted commas because they are using population models. I will use the nomenclature here for convenience to describe rate heterogeneity between branches, not within a branch (last example).

A "local clock" has mutation rate heterogeniety across the tree whilst a global clock is completely homogenious. You can't say how is a "local clock" related to global clock because there are numerous different models.

Global clocks never hold for biological data because there are always some taxa evolving at a faster rate than others. This is particularly true in infectious disease.

One historic model used in traditional local clocks (no inverted commas) was autocorrelation, where a local clock would relate to a global clock on a "clade by clade" basis. That is clades are correlated by their mutation rate, so some clades evolve faster and other clades slower. Statastically an autocorrelated clock always outperforms a global clock under a likelihood ratio test.

I remember autocorrelation clearly because I implemented a variant of this approach, so I needed to understand how they did it, rather than click a few buttons on a GUI, which is what happens in Beast.

Beast has at least 3 population models and off-hand I cannot remember the population scenarios. I do remember the population scenarios made sense and secondly they can give 10-fold differences in rate in the data set I looked at, You need to read on the background theory for each model and assess what is more biologically appropriate for a given context. One model (from dim memories) reminded me of a classic population growth equation, i.e. an exponential. Population growth and collapse in phylodynamics would be informative of the rate, but again you need to go to the underlying theory.

The molecular clock theory that the rate of evolution towards the tips of trees is much faster than at the deeper branches is the sort of model that could be implemented in a MCMC clock. In the latter case the relationship between a global clock and the rate heterogeneity clock (which is a better description in that model than "local clock") is quite different from global vs. "local clock".

The fundamental issue with clocks is exemplified in the 10-fold difference issue. The lesson that is important to me is very easy to get a "wrong answer". The lesson important to you is if two MCMC models can't agree the same rate how can both obtain the same answer as a global clock? I am personally a fan of maximum heterogeneity, but in this day and age they will ask you to perform a Beast model to support it.


HMM Profile from Convergence of retrotransposons in oomycetes and plants

The alignment they used to build the HMM is included as supplementary file 8 in the paper you cite:

  Multiple amino acid sequence alignment of RNH genes and domains from diverse taxa constructed and used for the phylogenetic reconstruction in the present study. (TXT 45 kb)

So, you can download that alignment and then use a tool like hmmbuild (part of HMMER) to build the HMM profile:

hmmbuild model.hmm 13100_2017_87_MOESM8_ESM.txt

The command above will create the file model.hmm which contains the HMM built from the alignment in the supplementary file 13100_2017_87_MOESM8_ESM.txt. I am including the first few lines here (this is not the entire file, don't use this as the model—build it yourself. I only include it to illustrate the format):

HMMER3/f [3.1b2 | February 2015]
NAME  13100_2017_87_MOESM8_ESM
LENG  121
ALPH  amino
RF    no
MM    no
CONS  yes
CS    no
MAP   yes
DATE  Mon Dec 11 13:33:29 2017
NSEQ  165
EFFN  5.697556
CKSUM 2893352190
STATS LOCAL MSV       -9.8975  0.71332
STATS LOCAL VITERBI  -10.6450  0.71332
STATS LOCAL FORWARD   -3.9247  0.71332
HMM          A        C        D        E        F        G        H        I        K        L        M        N        P        Q        R        S        T        V        W        Y
            m->m     m->i     m->d     i->m     i->i     d->m     d->d
  COMPO   2.34473  4.21240  2.85435  2.69203  3.54228  2.79153  3.54812  3.01140  2.72630  2.48370  3.82225  2.94269  3.49156  3.05882  2.84632  2.63343  2.86306  2.72633  4.27012  3.56181
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.13773  5.45955  2.08419  0.61958  0.77255  0.00000        *
      1   2.99488  4.57852  4.44659  3.29307  1.80597  4.16203  3.88788  2.87265  3.75733  2.38594  3.30344  3.33633  4.53284  3.64324  3.69630  2.90350  3.30742  2.83413  4.36991  1.50673      1 y - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00723  5.32905  6.05140  0.61958  0.77255  0.45275  1.01027
      2   2.66602  2.17750  5.01715  4.40788  1.72617  4.29647  4.62737  2.89579  4.21116  2.71668  3.64495  4.43195  4.66188  4.35829  3.99238  3.10363  2.02207  1.89528  5.14692  3.02155      2 f - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00661  5.41881  6.14115  0.61958  0.77255  0.61842  0.77391


Output of Seurat FindAllMarkers parameters

pct.1– The percentage of cells where the gene is detected in the first group
p_val_adj– Adjusted p-value, based on bonferroni correction using all genes in the dataset. This is also known as a false discovery rate (FDR) adjusted p-value. An adjusted p-value of 1.00 means that after correcting for multiple testing, there is a 100% chance that the result (the logFC here) is due to chance.
You are correct. The expression of Cttnbp2 donwn-regulated by that factor. If you were to look in the table at the cluster 1 DE genes you would see Cttnbp2 with a logFC value of 1.35264. When you are looking at genes differentailly expressed in only two clusters, the sign is the only thing that changes.

Vignette posted above by @heathobrien.


What Ensembl genome version should I use for alignments? (e.g. toplevel.fa vs. primary_assembly.fa)

There's rarely a good reason to use a hard-masked genome (sometimes for blast, but that's it). For that reason, we use soft-masked genomes, which only have the benefit of showing roughly where repeats are (we never make use of this for our *-seq experiments, but it's there in case we ever want to).

For primary vs. toplevel, very few aligners can properly handle additional haplotypes. If you happen to be using BWA, then the toplevel assembly would benefit you. For STAR/hisat2/bowtie2/BBmap/etc. the haplotypes will just cause you problems due to increasing multimapper rates incorrectly. Note that none of these actually use soft-masking.


How to load and split my two individual datasets when integrating datasets in Seurat?

Seurat is expecting individual datasets to be normalized separately prior to data integration. In that respect lists of objects corresponding to different datasets are handy to manipulate each object/dataset individually.

In your case you do not need to do any splitting, you have an object corresponding to each of your datasets.

If you would have started with a single count matrix of all datasets then you would need to split your data based on batches/experimental condition/...


Writing a python module to run Infernal on a list of fasta files

I think what you want to do is:

import subprocess
import os
import glob

genomes = glob.glob(os.path.join("/path/to/genomes/", "*.fa"))

for gfile in genomes:
    subprocess.Popen(args=["cmsearch", "/path/to/covariancemodel/", gfile], shell=False)

each pass through the for loop will take a one of the genome files and pass it into the cmsearch command.

Note:

1) glob.glob returns a list so no need to cast it again - https://docs.python.org/2/library/glob.html


What are the available cloud computing services for bioinformatics?

I have trialed the free version of InsideDNA, and these were my notes:

Cost: \$225/month for a team of 5 with 50TB storage or \$45/month with 10TB storage for individuals (assuming 6 month package: https://insidedna.me/pricing).
Software installed: Around 600 bioinformatic tools available and standard command line tools; some popular tools missing (like CD-HIT), but should be possible to install on request.
Jobs: Maximum of 32 CPUs and 208 RAM per job submission. Test jobs generally worked, although a larger test job failed.
Other points: Command line was sometimes slow, wget queries were slow, and scp was blocked. However, these may be resolvable issues.

Overall, I felt InsideDNA could be useful for groups without their own computational infrastructure and could be used for easily sharing resources between groups. The packages on offer seem not expensive, but I had a few issues, and I don't know how good their sys admin support would be.

I have not used the Amazon service, so can't comment beyond the details on their website. Also there are a few alternative companies, such as Genestack and DNAnexus, but I haven't directly tested them either.


Extracting concatenated metadata from a column in R

In general you want to:

Clean your text
Split the features by semicolons
Split keys and values by '='
Create a temporary data frame with keys and value
Reshape

library(tidyverse)
library(urltools)
library(textclean)
test_df <- data.frame(pos=c(3002983, 3002881, 3010946, 3021775), metadata=c("ID=SRX661585;Name=Ctcf%20(@%20Brain);Title=GSM1446329:%20CTCF%20brain%20ChIPSeq%3B%20Mus%20musculus%20x%20Mus%20spretus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Mouse%20adult%20hybrid%20brain;tissue=Adult%20whole%20brain;strain=BL6/spretus;chip%20antibody=CTCF,%20MillIpore%2007-729;", "ID=SRX2957352;Name=Smc1a%20(@%20AtT-20);Title=GSM2684855:%20SMC1%20ChIPseq%20in%20AtT-20%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=AtT-20%20cells;cell%20line=AtT-20%20cells;cell%20type=pituitary%20corticotroph%20cell;passages=10-20;chip%20antibody=Bethyl%20A300-055A;", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);"))

tbl <- apply(test_df, 1, function(row){
  # Split by smicolons and clean html tags
  n <- replace_html(unlist(strsplit(row[2], ";")))
  # Decode urlencoded characters, mostly space
  n <- url_decode(n)
  # Split into keys and values
  s <- do.call(rbind, str_split(n, "="))
  # Get keys and values into vectors and clean up leading and trailing whitespace
  v <- s[, 2]
  v <- str_replace_all(v, "(^ +)|( +$)", "")
  k <- s[, 1]
  k <- str_replace_all(k, "(^ +)|( +$)", "")
  # Create final data frame with original columns (pos), and also columns from
  # keys and values
  tibble(pos = row[1], key = k, value = v)
})

# Concatenate list into tibble
(tbl <- do.call(rbind, tbl))

This tibble looks like:
# A tibble: 37 x 3
   pos     key           value
   <chr>   <chr>         <chr>
 1 3002983 ID            SRX661585
 2 3002983 Name          Ctcf (@ Brain)
 3 3002983 Title         GSM1446329: CTCF brain ChIPSeq; Mus musculus x Mus spretus; ChIP-Seq
 4 3002983 Cell group    Neural
 5 3002983 source_name   Mouse adult hybrid brain
 6 3002983 tissue        Adult whole brain
 7 3002983 strain        BL6/spretus
 8 3002983 chip antibody CTCF, MillIpore 07-729
 9 3002881 ID            SRX2957352
10 3002881 Name          Smc1a (@ AtT-20)
# … with 27 more rows

Using this, you can:
# Pivot wider
tbl %>%
  pivot_wider(names_from = key, values_from = value)

Which nets you:
# A tibble: 4 x 15
  pos   ID    Name  Title `Cell group` source_name tissue strain `chip antibody` `cell line` `cell type`
  <chr> <chr> <chr> <chr> <chr>        <chr>       <chr>  <chr>  <chr>           <chr>       <chr>
1 3002… SRX6… Ctcf… GSM1… Neural       Mouse adul… Adult… BL6/s… CTCF, MillIpor… NA          NA
2 3002… SRX2… Smc1… GSM2… Neural       AtT-20 cel… NA     NA     Bethyl A300-05… AtT-20 cel… pituitary …
3 3010… SRX5… Mecp… GSM3… Neural       Olfactory … NA     NA     MeCP2 antibody… NA          Olfactory …
4 3021… SRX5… Mecp… GSM3… Neural       Olfactory … NA     NA     MeCP2 antibody… NA          Olfactory …
# … with 4 more variables: passages <chr>, `strain/background` <chr>, `genotype/variation` <chr>, age <chr>


how to search NCBI to find all available reviewed accession numbers for a protein family

A simple approach is to search for the gene name in NCBI Gene portal but this can miss a few cases and often include genes that are unrelated because the gene symbol is used by other unrelated genes as well.

For a cleaner dataset, you can use the NCBI Orthologs as follows. Search for human ALDP in the NCBI homepage with All Databases selected in the dropdown menu as shown below:

After you click the 'Search' button, you will be directed to a result page that looks like this:

Click on the 'Orthologs' button and you will be directed to a page with the list of proteins that are similar to the one you have started from. Note, these may not be orthologs in the truest sense; see this description about how they are calculated.

You can search for the taxonomic group that you are interested in, pick and choose the organisms of your interest and download the table in CSV format. It will have the NCBI GeneID, RNA and protein accessions.


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

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

Using wget to download the accession used as example:

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


Convert a BAM file from one reference to another?

You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but since no one has added NCBI names (yet, maybe I'll add them now) you're currently out of luck there.

Anyway, the quickest way to do what you want is to samtools view -H foo.bam > header to get the BAM header and then change each NCBI "chromosome name" to its corresponding UCSC chromosome name. DO NOT REORDER THE LINES! You can then use samtools reheader and be done.

Why, you might ask, would this work? The answer is that chromosome/contig names in BAM files aren't stored in each alignment. Rather, the names are stored in a list in the header and each alignment just contains the integer index into that list (read group IDs are similar, for what it's worth). This also leads to the warning above against reordering entries, since that's a VERY convenient way to start swapping alignments between chromosomes.

As an aside, you'd be well served switching to Gencode or Ensembl chromosome names, they're rather more coherent than the something_random mess that's present in hg19 from UCSC.

Update: Because I'm nice, here is the conversion between NCBI and UCSC. Note that if you have any alignments to patches that there is simply no UCSC equivalent. One of the many reasons not to use UCSC (avoid their annotations too).


You can't (easily) use a single design, since clone is nested within diabetic status, so you were correct using two separate designs. Your designs are correct.

If you really wanted to use a single design for some reason, it'd be ~clone+structure and for the latter comparison you'd end up using a rather complicated contrast (something like c(rep(c(1,-1), 3), 0)). I really can't recommend that unless you have a very compelling reason, since it's vastly easier to think in terms of factors for your experiment and you're not hurting anything by fitting twice.


How to find which scales are available in SeqUtil of Biopython?

Use this line of code:

print(analysed_seq.protein_scale(analysed_seq.get_amino_acids_percent(), window=9, edge=0.4))


If you want to call heterozygous SNPs and INDELs, you may try polyphred. I believe there are a couple of other similar tools, but I have forgotten their names.

  How can evidence of (het) indels like this be distinguished from noise in the chromatogram?

Your two points are correct. Basically, at some position, you start to see two streaks of overlapping peaks. Once identified, the two streaks can be aligned to each other if you shift one of them. This is a strong constraint, but it is hard to tell by eye because it is not obvious how to assign peaks to streaks without running an automatic tool. Doing this is even more difficult given noisy traces like your last figure. I think you'd better run polyphred or similar.

  the putative indel seems to start in a homopolymer region—cast doubt on this?

The INDEL mutation rate in long homopolymer runs is orders of magnitude higher than elsewhere in the genome. This has two implications: 1) you expect to see more homopolymer indels; 2) PCR polymerase is more likely to introduce errors in homopolymers––a higher fraction of homopolymer indels tend to be wrong.


How to create artificial or dummy phenotype data

Using wakefield package it is pretty easy, here is an example from GitHub page. This will create 500 dummy variables for id, race, age, etc. :

library(wakefield)

myDummyPheno <- r_data_frame(
  n = 500,
  id,
  race,
  age,
  sex,
  hour,
  iq,
  height,
  died
)

# result
myDummyPheno

# # A tibble: 500 x 8
#    ID    Race        Age Sex    Hour           IQ Height Died
#   <chr>  <fct>     <int> <fct>  <S3: times> <dbl>  <dbl> <lgl>
#  1 001   Bi-Racial    27 Male   00:00:00     93.0   71.0 F
#  2 002   White        28 Female 00:00:00    104     64.0 T
#  3 003   Black        28 Male   00:00:00    111     67.0 T
#  4 004   Hispanic     35 Male   00:00:00     90.0   67.0 T
#  5 005   Hispanic     21 Female 00:00:00     94.0   70.0 T
#  6 006   Black        35 Male   00:00:00     74.0   73.0 F
#  7 007   Black        27 Male   00:00:00    102     68.0 T
#  8 008   White        28 Female 00:00:00     95.0   70.0 F
#  9 009   White        30 Female 00:00:00     80.0   71.0 F
# 10 010   White        30 Female 00:00:00    114     66.0 T


Software to produce a table of post-translational modifications from a peptide list

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

Phosphorylation: GPS

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

Mitochondrial import: Mitoprot

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

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

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

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


Calculating Synonymous and non-synonymous

Woaa that is complicated data.

The reason it fails is easy, correcting it is harder. It fails because the open reading frame is lost during in silico translation, e.g. several gaps have 4 indels and the final nucleotides after the final indel are always 4. Thus it will invariably encounter a stop codon OR it will simply refuse to deal with an indel if it is incorporated into a triplet codon. There is no triplet codon model for an indel so the whole calculations stops.

Correcting data that complex ... the easiest way to do it is align you entire data set using amino acids, e.g. muscle. Have an identical nucleotide data set. MAKE SURE the names and positions are identical between your amino acid and nucleotide data, as a programmer this isn't hard to do however, everyone has their own scripts to do this. If you don't do this you next step will fall over.

Use a amino acid to nucleotide alignment translator. I use transalign (I think its EMBOSS). This will give you an inframe alignment. I noted many/all of you start codons are ATG which is a good sign. Check the nucleotide output and make sure there is no overhangs of the triplet codon at the terminal 3' end and check in Jalview or whatever your alignment editor of choice is that the translation is okay.

This will alignment work in the synonymous/non-synonymous analysis and all subsequent/complimentary analyses.

If you have an intron in protein your you'll need to remove it (hopefully not).

Finally you need to be careful of '3rd codon saturation' in your data because that heavily affects your trees and its certainly a risk given the heterogeneity of the indel structure, but first I would get an inframe alignment. There are checks for saturation (I'm sure even R can do this).


Error: FASTA file not found using building kallisto index

As mentioned in a comment, your working directory lacks the fasta file.


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.


It is important to train an error model as NanoSim does, as we do not fully understand the error processes involved in both the nanopore sequencing process and the basecalling process. Any sort of read simulator that does not use a species-specific model is just not going to produce realistic simulated reads.

As far as datasets go, I would recommend looking at this assembly and the associated data uploaded on EMBL:

  Michael, Todd P., et al. "High contiguity Arabidopsis thaliana
  genome assembly with a single nanopore flow cell." Nature
  communications 9.1 (2018): 541.

There are several data types there, as well as an assembly to train the models on.

Is it OK to train a model using reads from one individual and a reference assembly from another individual?

Generally, it is not good to map reads from one individual to the reference genome of another individual since this introduces something called reference bias. When we align reads to a reference genome generated from the same reads, or just from other reads from the same individual we remove this problem. All true reads will map to correctly assembled parts of the genome.

Why can we train a model for nanopore reads on a reference generated from ... nanopore reads? Won't we get a bad model for indels?

Yes, well, maybe. If the reference genome has been corrected to remove indels using a tool like pilon, then you can generally presume that the reference genome has very few to no indels present as a result of sequencing error from the Nanopore/PacBio reads. As a result, the model derived from aligning reads should have the correct error profile for indels.

Again, it is important that the long-read-assembly was corrected with Illumina shotgun reads from the same individual from which the long reads were derived. Otherwise the corrected reference would have a chimaera of indels and SNPs from multiple individuals.


Maybe comm, which is part of GNU coreutils, is more efficient.

comm writes to standard output lines that are common, and lines that are unique, to two input files; a file name of ‘-’ means standard input. Synopsis:
comm [option]… file1 file2

In your case you could sort both files based on the first 3 columns:
file1=file1.vcf
file2=file2.vcf

comm -12 <(awk '{print $1,$2,$3}' $file1 | sort) \
         <(awk '{print $1,$2,$3}' $file2 | sort) > matches

# result:
cat matches
1 10366 rs58108140
1 51954 rs185832753
#CHROM POS ID

This could be of course polished by removing the header.
Also if you want the original lines as output, you can grep the matches to the original files:
grep -f matches file1
#CHROM POS ID REF ALT QUAL
1 10366 rs58108140 G A
1 51954 rs185832753 G C

Edit: In the first version, intermediate temp files were used, but as suggested by terdon in the comments, this can be avoided using Process Substitution allowing the <(..) contstruct.


How to write a hash function for canonical kmers?

You can convert any string hash function to a "canonical" DNA string hash function.

Given a DNA string $s$, let $\overline{s}$ be its Watson-Crick reverse complement. Suppose $h:\Sigma^*\to\mathbb{Z}$ is an arbitrary string hash function. You can define
$$
\tilde{h}(s)\triangleq \min\{h(s),h(\overline{s})\}
$$
Then for any DNA string $s$
$$
\tilde{h}(s)=\tilde{h}(\overline{s})=\tilde{h}({\rm canonical}(s|h))
$$

PS: There are more than one ways to define "canonical" hash function. For example, let $h_0(\cdot)$ be the 2-bit encoding hash function. Given any integer hash function $g(\cdot)$, we can define:
$$
h'(s)\triangleq g\big[\min\{h_0(s),h_0(\overline{s})\}\big]
$$
Then we still have
$$
h'(s)=h'(\overline{s})
$$
You can also replace $\min\{\}$ with XOR or plus to get new definitions.


Interpretation of the violin plots from sc-RNA-seq

The black dots represent the values for individual cells. The "nGene" plot (the first one) shows the number of detected genes for every cell.

The red shape shows the distribution of the data. For the "nGene" plot, you can see that the average number of genes per cell is about 900 and most of the cells have roughly around 700-1100 genes.


Merging / Concatentating BLAST Databases

There's an app for that, apparently.

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

  Aggregate existing BLAST databases

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

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

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


Confirmation about PERMANOVA and Variance partitioning (permutest)

Permutational Multivariate Analysis of Variance (PERMANOVA) assumes no distribution which means it is a non-parameteric statistical test. According to the assumption of PERMANOVA test, each group must have same dispersion in the data. Before performing any statistical test, we must check the underlying assumptions. Variance partitioning (permutest) is a confirmatory test to check the assumption (similar dispersion). You can use permutest.betadisper for permutest. After that you can perform PERMANOVA test using adonis().
Interpretation of PERMANOVA test –

If p value is small (<0.05 ) we reject the null hypothesis (No difference between the groups) otherwise we accept the null hypothesis (Significant difference between the groups).

Coefficient of determination, R-squared value indicates the variation in the distances being explained by groups under the study. Higher value indicates a better fit.

You can further look into the confidence interval for each group.
Hope this confirms the answer to your question.


Why do reverse radtags have different start points in radtag sequencing?

Assume that after digestion you're left with the following sequence:

5' AGCTAGCTAGCTAGCTGACGTCGTAC 3'
3' GTACGACGTCAGCTAGCTAGCTAGCT 5'

Then one read would start with AGCTAG and the other would start with TCGATC.

The barcodes (what are used for multiplexing) are outside of this and part of the illumina adapters. They're sequences in completely different reads. This is nicely illustrated below:

Index 1 and 2 don't actually need to be present. Even if they are, whether the sequencer actually sequences them is a user setting.

For those curious, the order of the reads as they're sequenced are:

read #1
Index 1 (optional)
Index 2 (optional)
read #2 (optional)

If you ever get more than 2 reads from a sequencing provider, then it's likely that they included the index sequence (I have no clue why, but some sequencing centers like to do this...even though it confuses the heck out of their users).


Why does the SARS-Cov2 coronavirus genome end in aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa (33 a's)?

Good observation! The 3' poly(A) tail is actually a very common feature of positive-strand RNA viruses, including coronaviruses and picornaviruses.
For coronaviruses in particular, we know that the poly(A) tail is required for replication, functioning in conjunction with the 3' untranslated region (UTR) as a cis-acting signal for negative strand synthesis and attachment to the ribosome during translation. Mutants lacking the poly(A) tail are severely compromised in replication. Jeannie Spagnolo and Brenda Hogue report:

The 3′ poly (A) tail plays an important, but as yet undefined role in Coronavirus genome replication. To further examine the requirement for the Coronavirus poly(A) tail, we created truncated poly(A) mutant defective interfering (DI) RNAs and observed the effects on replication. Bovine Coronavirus (BCV) and mouse hepatitis Coronavirus A59 (MHV-A59) DI RNAs with tails of 5 or 10 A residues were replicated, albeit at delayed kinetics as compared to DI RNAs with wild type tail lengths (>50 A residues). A BCV DI RNA lacking a poly(A) tail was unable to replicate; however, a MHV DI lacking a tail did replicate following multiple virus passages. Poly(A) tail extension/repair was concurrent with robust replication of the tail mutants. Binding of the host factor poly(A)- binding protein (PABP) appeared to correlate with the ability of DI RNAs to be replicated. Poly(A) tail mutants that were compromised for replication, or that were unable to replicate at all exhibited less in vitro PABP interaction. The data support the importance of the poly(A) tail in Coronavirus replication and further delineate the minimal requirements for viral genome propagation.
Spagnolo J.F., Hogue B.G. (2001) Requirement of the Poly(A) Tail in Coronavirus Genome Replication. In: Lavi E., Weiss S.R., Hingley S.T. (eds) The Nidoviruses. Advances in Experimental Medicine and Biology, vol 494. Springer, Boston, MA

Yu-Hui Peng et al. also report that the length of the poly(A) tail is regulated during infection:

Similar to eukaryotic mRNA, the positive-strand coronavirus genome of ~30 kilobases is 5’-capped and 3’-polyadenylated. It has been demonstrated that the length of the coronaviral poly(A) tail is not static but regulated during infection; however, little is known regarding the factors involved in coronaviral polyadenylation and its regulation. Here, we show that during infection, the level of coronavirus poly(A) tail lengthening depends on the initial length upon infection and that the minimum length to initiate lengthening may lie between 5 and 9 nucleotides. By mutagenesis analysis, it was found that (i) the hexamer AGUAAA and poly(A) tail are two important elements responsible for synthesis of the coronavirus poly(A) tail and may function in concert to accomplish polyadenylation and (ii) the function of the hexamer AGUAAA in coronaviral polyadenylation is position dependent. Based on these findings, we propose a process for how the coronaviral poly(A) tail is synthesized and undergoes variation. Our results provide the first genetic evidence to gain insight into coronaviral polyadenylation.
Peng Y-H, Lin C-H, Lin C-N, Lo C-Y, Tsai T-L, Wu H-Y (2016) Characterization of the Role of Hexamer AGUAAA and Poly(A) Tail in Coronavirus Polyadenylation. PLoS ONE 11(10): e0165077

This builds upon prior work by Hung-Yi Wu et al, which showed that the coronaviral 3' poly(A) tail is approximately 65 nucleotides in length in both genomic and sgmRNAs at peak viral RNA synthesis, and also observed that the precise length varied throughout infection. Most interestingly, they report:

Functional analyses of poly(A) tail length on specific viral RNA species, furthermore, revealed that translation, in vivo, of RNAs with the longer poly(A) tail was enhanced over those with the shorter poly(A). Although the mechanisms by which the tail lengths vary is unknown, experimental results together suggest that the length of the poly(A) and poly(U) tails is regulated. One potential function of regulated poly(A) tail length might be that for the coronavirus genome a longer poly(A) favors translation. The regulation of coronavirus translation by poly(A) tail length resembles that during embryonal development suggesting there may be mechanistic parallels.
Wu HY, Ke TY, Liao WY, Chang NY. Regulation of coronaviral poly(A) tail length during infection. PLoS One. 2013;8(7):e70548. Published 2013 Jul 29. doi:10.1371/journal.pone.0070548

It's also worth pointing out that poly(A) tails at the 3' end of RNA are not an unusual feature of viruses. Eukaryotic mRNA almost always contains poly(A) tails, which are added post-transcriptionally in a process known as polyadenylation. It should not therefore be surprising that positive-strand RNA viruses would have poly(A) tails as well. In eukaryotic mRNA, the central sequence motif for identifying a polyadenylation region is AAUAAA, identified way back in the 1970s, with more recent research confirming its ubiquity. Proudfoot 2011 is a nice review article on poly(A) signals in eukaryotic mRNA.


How can I systematically detect unknown barcode/adapter sequences within a set of samples?

You mention that FastQC "fails to find the actual adapter sequences" - I guess you mean in the Adapter Sequence Contamination plot. However, the kmer and Sequence Content Plots are often useful even when the former fails. I've used these in the past - you can sometimes just read off the adapter sequence from the start of the Sequence Content Plot (or at least see how many bases to trim).


How to map reads shorter than 32bp with minimap2?

Minimap2 doesn't work well with 32bp reads by design. Older mappers like bowtie1 and bwa-aln will be better.


Adding to what benn pointed out, at this line:

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)

You specified that the counts, can be a explained by Genotype (being Control or Viral) or Type (RNAseq or Ribo) or an interaction effect

What you need is this coefficient, Genotype:Type. To use it, go

res <- results(dds,name="GenotypeViral.TypeRNA")

Based on your setup, you are always comparing RNA to Ribo (see "Type_RNA_vs_Ribo" in resultsNames). If changes to the abundance of a mRNA is similar in RNAseq or Ribo, then the coefficient of the interaction term will be zero and it will be zero in the log2foldchange column in results(res). A positive fold change indicates the log2foldchange between RNA and Ribo is higher in Viral.

If you would like to compare ribo versus and rna-seq, then set RNA as the reference by:

coldata$Type = factor(coldata$Type,levels=c("RNA","Ribo"))
# run deseq again
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)
dds <- DESeq(dds)

With this you will compare by default Ribo versus RNA, and the interaction term will tell you how much more enriched or depleted the ribo fraction is, in the Viral samples.

Would recommend checking the vignette out here:

http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts


I would recommend to do this in R, using biomaRt.

For example if you want all the SNPs from IDH1 (ENSG00000138413) gene:

library(biomaRt)

ensembl<- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl",mart=ensembl)

IDH1 <- getBM(attributes=c('ensembl_gene_id','variation_name'), filters = 'ensembl_gene_id', values ="ENSG00000138413", mart = ensembl)

IDH1

You can add much more information, take a look at the help page of biomaRt.


its easy replace match with sub but please stop using list as a variable ... mylist is fine.
This might work
   mylist = re.sub(r'H149xcV_\w+_\w+_\w+', gene_id)

otherwise
myregex = re.compile('_\w+\s+.*')
fastaid = myregex.sub('', myfile)

OR from @MaximilianPress
myregex2 = re.compile('_\w+\n') # or myregex2 = re.compile('_\w+$')
fastaid2 = myregex.sub('\n', myfile) # or fasaid2 = myregex.sub('', myfile)

The above will work ... As with all my code I never proof it ..


The .main and .sample file files usually come with default installations and may have been removed by whoever setup this Galaxy instance. You can look at the originals online and just base your own file on those: https://github.com/galaxyproject/galaxy/blob/dev/lib/galaxy/config/sample/tool_conf.xml.sample


Where can I find summary data for how common certain mutation *types* are?

This dbsnp github repo seems to offer that information.

https://github.com/ncbi/dbsnp/tree/master/Build%20Announcements/151


How can I get plant Genome size, total number of genes, EST, and chromosomes number

Ensembl plants has a REST API, you might be able to get all these details (except possibly EST numbers) there: http://ensemblgenomes.org/info/access/rest

For example this query gives you the Arabidopsis genome size: http://rest.ensemblgenomes.org/info/genomes/arabidopsis_thaliana?content-type=application/json


I have a list of protein targets. How to get structural data (atomic coordinate) from PDB?

There is not a simple way and there are factor that you also need to consider before starting.

Are you using gene names or gene accession numbers/id (NCBI/Uniprot)?
If not, are the names the current standard ones (e.g. scraped names from old papers is a definite no)
PDB residue index is rarely the same as protein index.
Do you want a perfect match or do you want say 80% homology?
Which PDBs do you want for a protein?

Personal example

For a website I have developed I have a route that gets all the PDBs for a given gene. (enter link description here). The way I did it is that I parsed Uniprot using a custom parser and made a dictionary/DB for each species with many synonyms for a gene going to a Uniprot ID and then from there I have another dictionary/DB that goes from Uniprot to PDBs if any. I went for Uniprot as I needed other data too (e.g. features) to present to the user. Including a renumbering option.
All this is way overkill for what you need. However, if you specify the species I can easily share a JSON file online.

SIFTs

If point 3 is important and you have Uniprot IDs I would suggest going with SIFTs. It has the offsets and has mapping for Uniprot IDs (pdb_chain_uniprot.tsv).

Gene names

In terms of mapping gene names to NCBI/Uniprot/Ensembl IDs, there are mapping tables everywhere. But if you may have obscure synonyms, using the NCBI API is a good option.
Uniprot data is better for protein and has PDB codes and chains, but the API gives a XML reply which is a bit annoying (I wrote a parser for it, but I think there was going to be one in Biopython).

Homologues

If you would like homologues, say your organism is human and you are okay with mouse, you can either use the threaded SwissModel models or use the blast query in RCSB PDB (not the PDBe, which has better data normally) or NCBI set to PDB database.
Swissmodel actually contains also PDB structures, so that may be a good starting point if your organism is one of those where precompiled exists.

Parts

Some proteins, especially mammalian are formed by separate domains linked together, especially the docking scaffolds. In which case you may care for a specific model within the gene, hence the indexing issue and the potential getting features.
If you want a single PDB structure for each gene and do not care about the technical I have some Pymol-using py3 code to fuse models together with the N and C termini aligned on line. But I strongly discourage that and would argue for protein feature based routes.


How to get log2 fold change of RNA-Seq data for time series experiment?

If you just want the change with time then use time as a continuous variable. The log2FC will then be the change per hour (you can use the default Wald test to extract it).

For time-point specific comparisons you'll instead need to treat the time points as groups and use contrasts as appropriate.


"Feature Barcoding" refers to using antibodies linked to specific oligos. The antibodies can then bind to cell surface markers of interest and the oligos they're conjugated with turned into normal scRNA-seq reads. This effectively means you can get a FACs profile for each cell without needing to run it through a FACs sorter and individually label it.

The "feature barcoding" nomenclature comes from 10X I think, but there are a few similar protocols (maybe they're exactly the same in the end, I'm not sure about the details) such as CITE-seq and REAP-seq. You can see a talk from 10X on this here.


A reliable fetcher of short read using SRA/ENA accession

Apparently there is no 100% reliable solution, but in the most cases, you should be able to use the enaBrowserTools python scripts (enaDataGet) to quite reliably download data.

The run ERR2135453 is a special case because the sample submission (ERS1939590) for it has been cancelled. The solution for these reads is to download the data manually from the ftp location.

This answer is based on kind assistance of Suran from ebi helpdesk.


Why my bim file doesn't match to my ped file as the Plink documentation suggests?

Responded in the plink2-users group (https://groups.google.com/forum/#!topic/plink2-users/zgJxdXxvdLo ).


How do I add a colour annotation bar to the heatmap generated by DoHeatmap function of Seurat v2?

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

From ?Seurat::DoHeatmap:

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

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

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

Description Condition
liver       Tumour
lung        Normal

Then you have to assign the colours

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

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

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


Confusion related to effective population size

That is a fantastic question and in my opinion a really important question in context. The Ne of Africa will be much bigger than the European population, however I'm less sure about the West African population which gave rise to the Afro-American population. I suspect it is still large because of the "out of Africa" theory.

My advice is to start by taking a complete look at the entire subject, from Tenesa et al (2007).  Ne is a really important parameter in any population study and extremely useful. I have feeling this could be complex, but I cannot hypotheses any further because I don't do human genetics (except for immunology).


Generate VCF from different .bam files with different chromosome names

The main issue here is that your bam files have different chromosome labels. (i.e. 1,2,3 vs. chr1,chr2,chr3) as you mentioned.

This suggests that the data was mapped to a slightly different version of the reference genome. As a first layer, you might want to double check that these reference chromosomes are cross compatible (e.g. both correspond to hg19 for instance). To find more information on the reference genome used, you can check the header information using "samtools view -H "

Assuming that the reference genomes are cross compatible, you can try converting it to sam, adding the missing "chr" label in a pipe, and converting it back to a bam file. However, this will also require modification of the header file which can be challenging in itself.

The easiest (and more rigorous) approach is probably just to convert the batch of bam files with the wrong chromosome labels back into fastq files. You will then want to remap them using the exact same reference genome as the other batch of bam files. After which, you should be able to perform variant calling on the combined batch of bam files.


featureCounts segmentation fault on Arch Linux

I resolved this problem by installing from source; I assume the precompiled binary was not compatible with the Arch Linux kernel.


How to download the whole BLAST nt database into a specific folder?

-O is to specify the name of the output file. You want the -P option.

So, from your example

wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz" -P output/


Calculating Z-score from logCPM values using edgeR

Yes this is nothing but a z-score. The scale function calculates column z-scores, but since you use t (transpose), it is actually performed on rows. There are many other ways to calculate z-scores, but if done correctly they will all give the same results.


Hmmmm .... you don't actually need pop ID to work for Structure, providing their are no logical errors (you get a wrong answer). The bit that must work is the a priori number of populations in the data set. You are then going to manually (or via a script) feed a range of a priori values into Structure, run each and plot the likelihood. Where the likelihood plateaus thats the number of populations contained within your data. This is the key output of Structure.
Personally, I would perform a Structure run on the test data and remove the pop ID column and check the answers (likelihood -lnL) are identical. Thereafter use an approach such as $F_{ST}$ to interrogate the exact members of each population delineation.

I think your test will work .... The only thing I can think of is the header comprises a space, This format just needs checking against the test data.


Are there databases or resources that contain expression of proteins or genes before and after cancer treatment?

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

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


Error while Importing fastq sequences to jMHC software

Yes, use .fasta , you can try to convert your fastq files with NextGen Workbench for example.


What command/invocation is used to generate NCBI 16SMicrobial blastdb

Running blastdbcmd -db foo -info provides a little information but I haven't seen anything which will report exactly how a blastdb was created.
A blastdb readme suggests that only the parse_seqids option has been added to the standard parameters.

For those from NCBI, the following makeblastdb commands are recommended:

For nucleotide fasta file:   makeblastdb -in input_db -dbtype nucl -parse_seqids
For protein fasta file:      makeblastdb -in input_db -dbtype prot -parse_seqids

The -parse_seqids requires NCBI fasta headers in the correct format (spec can be found here).