Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?

I'm trying to figure out why. This will be a longer read, Tl;dr at the end:

Doing a match

There is a good match of Q against LC074724 in 20170824:

$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
Q   LC074724.1  2581

Note that the bitscore of the LC074724 match is 2581.

But there is no match of LC074724 at all in 20170821:

$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724

The option -outfmt '6 qaccver saccver bitscore' tells BLAST to print the query accession number and version, the subject accession and version, and the bitscore of the match.

What are the best hits for 20170821

So what are the bitscores of the Q matches in 20170821?

$ blastn -outfmt '6 bitscore' -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | sort -nr | head
3301
2616
2590
2590
2563
2554
2547
2544
2540
2540

The bitscore of Q against LC074724 in 20170824 (2581) would have placed it in 5th position in the top 10 matches of Q, but LC074724 is not matched at all in 20170821! The lowest bitscores matched in 20170821 is 2354.

What's the difference between the databases?

Get the ids of the subjects in each database

$ blastdbcmd -entry all -db 20170824 | egrep '^>' | cut -c2- | sort > 20170824.ids
$ blastdbcmd -entry all -db 20170821 | egrep '^>' | cut -c2- | sort > 20170821.ids

$ wc -l 2017082?.ids
    9754 20170821.ids
    9111 20170824.ids
   18865 total

20170821 has 643 sequences that are not in 20170824:

$ comm -23 20170821.ids 20170824.ids | wc -l
     643

20170824 has no sequences that are not in 20170821:

$ comm -13 20170821.ids 20170824.ids | wc -l
       0

and the two have 9111 sequences in common:

$ comm -12 20170821.ids 20170824.ids | wc -l
    9111

These numbers match with the numbers from wc -l above.

So 20170821 is just 20170824 with an extra 643 sequences.

Maybe LC074724 matches, but with a low bit score?

BLAST only shows the top 500 matched subjects by default. Let's ask it to show more:

$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 -max_target_seqs 10000 > 20170821-10000-matches.txt

And....

$ head 20170821-10000-matches.txt
Q   Q   3301
Q   X234A   2616
Q   CS388973.1  2590
Q   DM059402.1  2590
Q   LC074724.1  2581
Q   KJ843188.1  2576
Q   AB697496.1  2576
Q   AB697499.1  2576
Q   AB697503.1  2576
Q   AB697508.1  2576

There it is! And it's in 5th position, just as it should be (based on the bitscores discussed above).

Just to confirm, taking away the -max_target_seqs option gets no match:

$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
$

Three confusingly named BLAST options

I've always found the names and descriptions of the following three BLAST command-line options very confusing.

From the Formatting options section of blastn -help:

  -num_descriptions =0>   Number of database sequences to show one-line descriptions for   Not applicable for outfmt > 4
  Default = `500'    * Incompatible with:  max_target_seqs

  -num_alignments =0>   Number of database sequences to show alignments for   Default = `250'    * Incompatible with:
  max_target_seqs

And from the Restrict search or results section:

  -max_target_seqs =1>   Maximum number of aligned sequences to keep   Not applicable for outfmt <= 4   Default = `500'    *
  Incompatible with:  num_descriptions, num_alignments

From the section names, it seems like only the latter option (-max_target_seqs) would have any effect on the search, and that the former two are just about what is displayed.

But... using either just -num_descriptions or -num_alignments with a big value also gets a match:

$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_descriptions 10000 | grep LC074724
LC074724.1  Hepatitis B virus DNA, complete genome, isolate: p621     2581    0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621

$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_alignments 10000 | grep LC074724
LC074724.1  Hepatitis B virus DNA, complete genome, isolate: p621     2581    0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621

So BLAST is apparently finding the LC074724 match even without the high -max_target_seqs option, but it's not displaying it unless there is a high -num_alignments or -num_descriptions value. (Keep in mind, this is the 5th best hit! Not somewhere >500)

Just to confirm, when none of the 3 options are given and no -outfmt option is either, the match is not found:

$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | grep LC074724
$

Tl;dr:

In summary, the options

 -evalue
 -max_target_seqs
 -num_alignments
 -num_descriptions

are NOT just filters applied to the results, but instead directly influence the results found.

This means you can End up with a bacterial sequence as BEST hit when using a value x, and with no bacterial sequences at all in your top20 when using value y.

To quote the NCBI user service:

  We don't consider this a bug, but I agree that we should document this
  possibility better. This can happen because limits, including max
  target sequences, are applied in an early ungapped phase of the
  algorithm, as well as later.

And

  The cutoffs are applied during an ungapped extension phase, and some
  alignments may improve in score during gapped extension. The safest
  approach is to relax (increase) the threshold and deal with a larger
  result set.

This was totally unexpected behaviour for me and is as far as I know is nowhere documented in the BLAST help or similar.


Significance of k-mer length in COVID-19 sequence analysis?

5bp produces a high resolution analysis, providing its sufficiently robust.

Its a sliding window of 5 bp wide and 1 bp step length which runs along the length of each genome. They are performing SNP analysis/classification. The original sequence they were examining was,

ATTAAAGGTTTAT

In any sliding window analysis the window size is arbitary towards the resolution and robustness of the analysis. The author believes that 5bp provides sufficiently robust information for a very high resolution classification/training analysis in their example.

The proof is the accuracy of the classification. I assume they have good results from 5bp. In ML the rationale is the general laws of ML (over/under training) and assessed on the test split in the data (probably 20%).

A key reason for the 1bp step is that this is the way ML works, you need a lot of data and by doing 5 bp sliding window at 1 bp step that will generate ~30 000 pieces of data for each sequence. So that's 3 million data points for 100 sequences.


Matching files by chromosome during wdl scatter gather

You could use the sub function to obtain the chromosome from the globbed file names and then use that to retrieve the correct LD matrices from a Map (see read_map or read_json).

eg.:

Map[String, String] ld_matrices = read_json(file_of_ld_matrices_by_chr)

scatter(gwas_for_chr in split.gwas_by_chr) {
    String chr = sub(sub(gwas_for_chr, prefix, ""), "\.assoc$", "")

    call sbayes.run {
        input:
            GCTB = gctb_executable_path,
            gwas = gwas_for_chr,
            ld_matrix = ld_matrices[chr],
            output_prefix = out
    }
}

The file_of_ld_matrices_by_chr would have to look like this:

{
  "chr1": "/path/to/file1",
  "chr2": "/path/to/file2"
}


For what bioinformatics tasks is Biopython more adapted than Bioperl?

Regading the perl vs python discussion, there is no final answer which language is better, but I have some advice for you:

Learn the language your colleagues or your advisor use. This way you are able to discuss your code with them and also get help if you run into problems.


Filter out other character after the genotype

Using sed:

sed 's,\([01][/|][01]\):[^ ]*,\1,g' file
#    ^  ^               ^     ^  ^
#    |  |               |     |  *--- applies to all found pattern
#    |  |               |     *------ print the content inside the parenthesis
#    |  |               *------------ any character except a space
#    |  *---------------------------- the pattern to be kept
#    *------------------------------- substitution

You also may want to format the output using column command:

sed 's,\([01][/|][01]\):[^ ]*,\1,g' file | column -t
#CHROM  POS     ID  REF   ALT   FATHER  MOTHER  DAUGHTER1  DAUGHTER2  DAUGHTER3   SON1  SON2
1       736689  .   T     A     0/0     1/0     0|0        0|1        0|1        0|0   0|1
1       758663  .   CGCC  TTAG  0/0     0/1     0|0        0|1        0|1        0|0   0|1

Note the comma separators in the sed command can replaced with other characters typically "/"


Where can I find the DNA sequence of the C. elegans rRNA cluster?

The info is actually in the gene_biotype column of my wormGenes data frame. I could select the rRNA genes using dplyr::filter(wormGenes, gene_biotype == 'rRNA'). Unfortunately, no description is available for the genes, so I will have to get it using the gene IDs.


Do people today with a disease caused by single gene mutation that is inheritable - all share a common ancestor?

That can happen, but it's more likely that the disease will be present in more than one family. Diseases with a single, small genetic cause will be more likely to flare up elsewhere (assuming the disease isn't too bad, and it wouldn't be if it's been in the family for 500 years).

However, what you will find is that in your family that disease is almost always present together with some other inherited genetic feature that is nearby on the genome (i.e. an inherited haplotype block). This association is less likely to be shared with other families.


Need help with Pairwise Alignment module in iterating over the alignment

aligner.score(ref_seq_1.seq , seq1.seq) doesn't produce an alignment. It just gives you the score of the alignment as float value.
Instead you need to run:
alignments = aligner.align(ref_seq_1.seq , seq1.seq)
Now alignments is an Bio.Align.PairwiseAlignments you can iterate over alignments or get the score with alignments.score which gives you the same as aligner.score(ref_seq_1.seq , seq1.seq)`
If you need more information, you could check the documentation.


The most intuitive explanation is  also explained in the background section:

  Conceptually, this methodology can be understood as a change in coordinate systems for gene expression data, from genes to gene sets.

Which I think it is explained in Figure 1 of the paper, which is also one of the most informative I found about the methods of an algorithm1.

I'll will use it to explain what's going on.

Fit a distribution function per gene

  For each gene expression profile ... a non-parametric kernel estimation ... function is performed.

I can't give more detail because I haven't understood it fully

Then they are normalized to make the ranks symmetric around zero

It is an inline formula between equations 2 and 3.
The Kolmogorov-Smirnov like random walk statistic is applied to those kernel estimations normalized distributions.

Which two distributions does it compare? The genes in a gene set and all the others (but I'm not totally sure)

The role of random walk is explained a bit later from your quote:

  [the random walk] produces a distribution over the genes to assess if the genes in the gene set are more likely to be found at either tail of the rank distribution

Which is later used to calculate the enrichment score by looking the maximum deviation from zero, or by summing the largest and the lowest deviations.

To address some other questions:

  GSVA calculates sample-wise gene set enrichment scores as a function of genes inside and outside the gene set, analogously to a competitive gene set test

1: Try finding which fittings and approximations are done inside limma


How can I export a full alignment from IGV as an image?

Running an IGV batch script helps increase the number of reads per image, compared with File | Save image, but it does not solve the issue by itself.

To run a script, first make a script file like so with 1 command per line:

snapshotDirectory ./
goto chr2:28582094-28582790
maxPanelHeight 10000
snapshot chr2_28582094_28582790.png

Run the script using Tools | Run batch script menu. The importat trick is the maxPanelHeight 10000 which is what makes this include more reads than the default snapshot.

You can also convert this into a little shell script:

#!/bin/bash
mkdir -p /tmp/igv
sleep 2
echo "snapshotDirectory /tmp/igv"
sleep 2
goto chr2:28582094-28582790
sleep 2
echo "maxPanelHeight 10000"
echo "snapshot chr2_28582094_28582790.png"

Then, just open IGV, load the relevant BAM file and run the script like this to connect to IGV via telnet (credit to Pierre Lindenbaum):

sh foo.sh |  telnet 127.0.0.1 60151


Tool for rna/lna melting temperature prediction

There are several programs for simulating/predicting DNA melting but you could try them for RNA.
You can look here http://bioinformatics.org/meltsim/wiki/Main/HomePage.


Alignment based vs reference-free (transcriptome analysis)?

First of all, I would emphasize that "alignment-free" quantification tools like Salmon and Kallisto are not reference-free. The basic difference between them and more traditional aligners is that they do not report a specific position (either in a genome or transcriptome) to which a read maps. However, their overall purpose is still to quantify the expression levels (or differences) of a known set of transcripts; hence, they require a reference (which could be arbitrarily defined).

The most important criterion for deciding which approach to use (and this is true of almost everything in genomics) is exactly what question you would like to answer. If you are primarily interested in quantifying and comparing expression of mature mRNA from known transcripts, then a transcriptome-based alignment may be fastest and best. However, you may miss potentially interesting features outside of those known transcripts, such as new isoforms, non-coding RNAs, or information about pre-mRNA levels, which can often be gleaned from intronic reads (see the EISA method).

This paper also has some good considerations about which tools may work best depending on the question you want to answer.

Finally, another fast and flexible aligner (which can be used with or without a reference transcriptome) is STAR.


Difference between de novo transcriptome assembly methods

Your intuition is correct. stringTie is just looking at clumps of alignments and how they might relate to each other (either due to spliced alignments or proximity). Trinity is doing the more computationally difficult task of finding parts of reads that overlap with each other and trying to link them together into longer sequences.

Whether it makes sense to use something like stringTie or not depends mostly on the quality of the reference sequence. If you're working on an organism with a high-quality reference sequence (pretty much anything with a lot of papers using it) then there's no reason to bother with Trinity. If your reference genome additionally has a high-quality annotation (human, mouse, etc.) then you should use it, as anything "novel" that the program is likely to output otherwise is probably just an artefact. If you either don't have an annotation or it's known to be highly problematic then it makes sense to omit it when using stringTie.

Without seeing evidence to the contrary, I don't buy the argument that using a good annotation is actually biasing the output in a problematic way. For well-studied organisms, anything reported as "novel" these days likely doesn't exist, or is otherwise so fleeting that it's biologically irrelevant. The arguments of bias sound like something that was come up with once someone didn't get as many "novel" results as they were hoping.


Let's simplify...
cmd.distance(xx,xxx) accepts two selector strings and is fine with and state 1. And best of all if you select multiple atoms it gives you the smallest distance!
So all you need to do is:
# measure
d_O = cmd.distance('UNL and element O', 'resi 312 and name ND2')
d_N = cmd.distance('UNL and element N', 'resi 113 and name OD2')
# write
f=open("dst4.txt", 'w')
f.write(f'Lig_N D113 {d_N:.2f}\n')
f.write(f'Lig_O N312 {d_O:.2f}\n')
f.close()

If you are using python 2 (if so do note pymol scripting is better in 3), then the second bit is
f=open("dst4.txt", 'w')
f.write('Lig_N D113 %.2f\n' % (d_N,))
f.write(f'Lig_O N312 %.2f\n' % (d_O,))
f.close()

Unsolicited coding tips
A few unsolicited tips... Your original script is rather over-complicated and obfuscated and uncommented, so for better legibility why not do:

cmd.iterate_state is always harder to follow than the for atom in get_model('xxxx').atoms synthax
unless you are stuck with python 2, use f-strings not the %s sprintf style formatting
the /x/x/x/x selector synthax is never as readable as the longer
resi/resn/name selection algebra.

Code golf is always tempting, but it can result in future-you getting stuck as happened here...


Some answers by the author Charles Grant via the MEME Suite Google group may be useful to you.

Via: https://groups.google.com/forum/#!topic/meme-suite/yNascbE8Tig/discussion

  Providing your own higher-order background model can greatly improve MEME's ability to discover motifs...

  The command-line utility fasta-get-markov, included in the MEME Suite download, is used to generate custom background Markov models. The input to fasta-get-markov is a FASTA file containing "background" sequences. Ideally, these "background" sequences will be different from the sequences you are analyzing with MEME, but as similar in nature as possible. For example, if you wanted to discover motifs in certain intergenic regions, you might use other sequence data from other intergenic regions to generate the background. The larger the set of "background" sequences is, the better the results will be...

  Typically, you should not specify an order larger than 3 for DNA sequences, or larger than 2 for protein sequences. However, if your input sequences contain higher-order non-random effects that are getting in the way of motif finding, you can follow the following "rules of thumb":

  Use a background model at least four orders less than the shortest motifs you are looking for. So, if you want to find motifs as short as six, I wouldn't use a model higher than order two.
  For an accurate model of order N, you need to use a FASTA file as input to fasta-get-markov with at least 10 times 4**(N+1) DNA characters in it.

  So,

  order-3 requires 2560 characters
  order 4 requires 10240 characters
  order 5 requires 40960 characters etc.

I'm not certain if that input requirement scales differently for proteins (i.e., 20^(N+1)). That may be worth following up on the Google group.

And via: https://groups.google.com/forum/#!topic/meme-suite/UiDhdY6WlK4/discussion

  You should choose the motif occurrence model based on your understanding of the experimental protocol that generated the sequence data and the underlying biology. MEME is trying to find short sequences that are statistically over-represented in your sequence data. To do this, it has to assume a model for how many occurrences of a motif there will be in each sequence. The nature of your experiment should be the basis for the model you choose. For example, data from a ChIP-Seq experiment will typically contain exactly one occurrence of a motif per sequence. On the other hand, sequences drawn from a database of upstream regions of co-regulated genes might contain any number of motif occurrences, including 0.

Via: https://groups.google.com/d/topic/meme-suite/QJiJsy1QxYk/discussion

  The MEME algorithm run time is cubic with respect to the number of input sequences, therefore, it is unsuitable for OOPS (only one per sequence) analyses that have more than 1000 sequences. For ZOOPS (zero or one per sequence) and ANR (any number of repeats) analyses, having more than 1000 sites are also intractable in terms of run time. For example, on a 3 GHz CPU, initialization of the p-value table takes about 10 min for 1000 sites/sequences in OOPS, about a day for 5000 sites, and about 1 week for 10,000 sites. Then the quadratic run time with respect to the total number of characters kicks in next.

Thus it seems like the ZOOPS and ANR models have similar runtime complexity, and use of the MPI-based parallel meme_p binary may be helpful, depending on your analysis:

  The parallel version of MEME scales up to about 128 processors. Please see http://www.sdsc.edu/~tbailey/papers/cabios96.pdf for a discussion of the parallel program. You must bear in mind however, that doubling the total number of sequences input to MEME means that you will need 8 times more processors for the job to finish in the same amount of time.


How to count the kmer occurrence in FASTA file considering overlapping and reverse complement?

There's a lot detail missing here, but at a guess could you try,
from Bio import SeqIO
from Bio.Seq import Seq
for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    reversec = record.seq.reverse_complement())
    revcount = reversec.count_overlap('AT')
    print (revcount + count_overlap('AT'))

Essentially count_overlap only counts the forward strand if this is the Biopython method. Not debugged the above.


What could be the reason for the samples not clustering?

Samples will only cluster by experimental group if the experimental effect is large enough that it's the primary source of variance between your samples. If that's not the case then you'll get results like you're seeing. Whether that's a problem or not ends up depending on how large of an effect you expect to see. If your groups are different stages in differentiation or another large effect, then you expect that samples will cluster into groups. If you're looking at something like some large unaffected tissue in a disease vs. normal state then seeing a cloud of points is quit reasonable.


Generating the reconstructed alignment from BAM

If all you want are the number of insertions/deletions, you can use this to extract that info:

from cigar import Cigar
c = c=Cigar('20M5I10M3I5M')
num_insertions = len(list(filter(lambda x: x[1] == 'I', c.items())))
insertion_length = sum([x[0] for x in list(filter(lambda x: x[1] == 'I', c.items()))])

Deletions would be the same thing with 'D' in place of 'I'


Changing the axis limits of ggplot objects

Axes for ggplots can be controlled via xlim() and ylim() functions. You can try:

your_ggplot_object +
  ylim(-5,20) # assuming that you want the gene labels to be readabe, not clear from your post


Can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence?

For DMPfold you might like to look at Figure 5B of the paper, which implies that performance starts to drop off after ~600 residues. Only proteins up to 500 residues were used in the training set, for example.

A protein of that length is almost certainly multi-domain. I would consider doing domain prediction; we have a tool DomPred, but many others are available. Then run the sequences of the domains individually and recombine them later.

I don't know about I-Tasser's ability to do this, perhaps see if they mention it in their papers.


How to obtain clusters of hierarchical heatmap when using Python?

Clustering like this is typically done with scipy. Here's the code we use in deepTools (original context here):

from scipy.cluster.hierarchy import fcluster, linkage
Z = linkage(some_matrix)  # You might want to set `method` and `metric`
groups = fcluster(Z, nGroups)  # You might want to set `criterion`

groups is then a vector containing the group assignment of the input matrix according to the clustering. As an aside, this is exactly the same as what's shown in the last link in your post, so I'm not sure how you didn't find it.

BTW, at least in python, hierarchical clustering is painfully slow with large datasets.


Can I get longer contigs by changing MEGAHIT settings?

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


protein secondary structure matching/comparison scoring

Hey i solve this simply by use of SOV alghoritm. https://www.researchgate.net/publication/319526145_A_modified_definition_of_Sov_a_segment-based_measure_for_protein_secondary_structure_prediction_assessment

cheers


split 1 column input into 5 column bed file

The following Python script seems to do the job.

#!/usr/bin/env python3

from __future__ import print_function
import re
import sys

for line in sys.stdin:
    if '_' in line:
        match = re.search(r'(\S+):(\d+)_(\d+)del([ACGT]+)ins([ACGT]+)', line)
        print(match.group(1), match.group(2), match.group(3), match.group(5), match.group(4), sep='\t')
    elif 'dup' in line or 'del' in line:
        match = re.search(r'(\S+):(\d+)(del|dup)([ACGT]+)', line)
        if match.group(3) == 'del':
            print(match.group(1), match.group(2), match.group(2), match.group(4), '-', sep='\t')
        else:
            print(match.group(1), match.group(2), match.group(2), '-', match.group(4), sep='\t')
    else:
        match = re.search(r'(\S+):(\d+)([ACGT]+)>([ACGT]+)', line)
        coord = int(match.group(2))
        print(match.group(1), coord - 1, coord, match.group(3), match.group(4), sep='\t')

Invoked like so on the command line.

[standage@lappy ~] $ ./transform < f1
chr7 140453144 140453145 A T
chr7 140453135 140453136 A T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -


We use kinship2 R package to visualize pedigree files (.fam plink files). For some reasons, some characters from plink ped aren't supported (like "-9" for unknown genre or affected status).
Here is the R script we use:

args<-commandArgs(TRUE)
pedPath<-args[1]
pdfFile<-args[2]

library(kinship2)
pedfile <- read.table(pedPath,sep="\t")
pedfile$V5 <- apply(pedfile,1,function(data){if(data[5] == 1 || data[5] ==2){data[5]}else{3}})
pedfile$affected <- apply(pedfile,1,function(data){if(as.numeric(data[6]) == 1){res <- data[6]}else{if(as.numeric(data[6]) == 2){res <- data[6]}else{res <- NA}}})

familyIDs <- as.character(unique(pedfile$V1))

pdf(file=pdfFile)
for(famID in familyIDs){
    ped <- pedigree(id=as.character(pedfile$V2[pedfile$V1==famID]),dadid=as.character(pedfile$V3[pedfile$V1==famID]),momid=as.character(pedfile$V4[pedfile$V1==famID]),sex=as.numeric(pedfile$V5[pedfile$V1==famID]),affected=as.numeric(pedfile$affected[pedfile$V1==famID]))
    plot.pedigree(ped,cex=0.7,align=TRUE)
    title(famID)
}

And we call it this way:

Rscript <kinshipScript.R> <plink ped file> <pdf output>


cellranger mkfastq with full path to --id flag

You should be able to use the --output-dir argument to specify a folder of your choice and use --id=demux_10x as before.


A Poisson or binomial would suffice.

I'm not at all sure the difference will significantly deviate from randomness.

Perl: 95% confidence interval of a non-parametric bootstrap for e.g. 1000 replicates is an alternative approach. I'd be more comfortable coding that solution in Perl than I would in Python, so is unlikely to be of much use here. Anyway, I would use a 2-dimensional array for your total data-set, the second dimension being coding/non-coding, I would then use the inbuilt random-number generator and multiply by the size of the array to randomly pluck out each element. 95% CI of the mean works by ranking the data and removing the top and bottom 2.5% and taking the interval for a single replicate (which is repeated 999 times more).

There will be a deep learning approach to this data, which would be easier to implement in Python, but I'd need to think about this because it is not simple classification.

Again my warning is that given the sample size, the result at best might be borderline. There are likely to be lots of interesting in your data, but coding/non-coding may not be one of them.


Which are the use cases for the methods for DE in Seurat

You can take a look at the recently published article: Bias, robustness and scalability in single-cell differential expression analysis.

  We evaluated 36 approaches using experimental and synthetic data and found considerable differences in the number and characteristics of the genes that are called differentially expressed. Prefiltering of lowly expressed genes has important effects, particularly for some of the methods developed for bulk RNA-seq data analysis. However, we found that bulk RNA-seq analysis methods do not generally perform worse than those developed specifically for scRNA-seq.

Then you can decide which method fits better to your type of scRNA-seq protocol. Since Seurat offers DESeq2 could be a good idea to use it, but this depends a lot on your data.


BERT Language Model and Gene Sequences - How Do I Relate Clusters of Sequences?

The UMAP results look very pretty :)

If I were to try to make sense of this blob, I'd try to do the following.

Pull out Gene Ontology (GO) terms for a list of candidate genes (you can try to subset Uniprot / Swissprot for some examples). Then I feed through some genes through your fragmented BERT kmer approach to see where they land. I'd be particularly interested in seeing where transmembrane proteins and enzymes land here.
I'd try to color this by introns vs exons (you should be able to extract this sort of information from your genome files). There is a good chance that the difference between introns / exons could explain some of these differences.
It maybe worth looking into genes that only encode RNAs vs genes that encode proteins.


Using data mining of papers in order to derive genomic connections

There are many databases that have used publication scraping for oncogenic gene fusions. There are publications for the individual methods they used for their scraping and aggregation.

COSMIC - http://cancer.sanger.ac.uk/cosmic
TICdb - http://www.unav.es/genetica/TICdb/
ChimerDB - http://ercsb.ewha.ac.kr/fusiongene/ -
https://www.ncbi.nlm.nih.gov/pubmed/27899563


How to search for high coverage SRA entries

If processing a huge table with lots of "false positives" is OK, then you could enter your search term, from top right choose "Send to" > "File" Format: "RunInfo".


Advantages of paired-end sequencing compared to single end

In paired-end sequencing, you have an expected fragment length. So if Read 1 maps to two different places, whereas Read 2 only maps to one place. Then you can determine which one is more likely for Read 1 to map to. It cannot be too far away from Read 2 given the expected fragment length.
It is not necessarily true that paired-end sequencing gets twice the amount of data. The total amount of data can be determined when running sequencing. If single-end has 20M reads, I can just stop sequencing at 20M as well for paired-end.


I think you can come close to this with cutadapt:

cutadapt \
  -a ${FWD_ADAPTER} -A ${REV_ADAPTER} \         ## adapter
  --max-n 0.1 \                                 ## do not allow > 10% Ns
  -q 5 \                                        ## remove bases with basequal < 5
  --minimum-length (readlength/2) \             ## minimum length after trimming
  -o read1_trimmed.fq -p read2_trimmed.fq \     ## outputs
  read1.fq read2.fq                             ## inputs

This will remove reads with more than 10% of bases being N and remove bases with base quality below 5. If you then set --minimum-length to half the read length, this should also remove reads where more than half of the bases had a base quality below 5 and therefore got trimmed away. Poor base quality in Illumina reads typically happens towards the 3' end so I think this little workaround is a sufficient proxy to achieve your goal. A base quality of 5 is very low in any case, so if sequencing went normal the typical quality should be notably higher. I would be surprised if this third filter is really necessary and will eliminate a notable number of reads.


Viral genome assembly using broad viral ngs pipeline?

Since you are already using the Broad Tools sets you can use Picard FastqToSam to make the conversion

As far a clipDb I am unfamiliar with that and a quick google search and look at the trimmomatic manual were unhelpful


How to specify resources for cluster in snakemake

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

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

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

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

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

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

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

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

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

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


Pandas automatically rounds GWAS P-value

The solution you likely want is here,

pd.set_eng_float_format(accuracy=x, use_eng_prefix=True)

x = whatever is required
The function set_eng_float_format has been moved around a bit and is now a top level function

You might be dealing with maximum likelihood, so convert to log likelihood -200 that sort of thing

Best idea

gwas_rs_id['logged'] = np.log(gwas_rs_id.marker)

First idea,

import decimal as D
....
gwas_rs_id['marker'] = gwas_rs_id['marker'].astype(D.decimal)

Last attempt,

pd.set_option('display.precision', x) # x whatever dp you need, this might work in combination with a log likelihood

Try a gwas_rs_id.dtype() check whether you're using float32, float64 or still in object. It should convert to scientific notion (E values) automatically.


This is possibly more of a molecular genetics question appropriate to biology.stackexchange.com, also note that it's usually better on stackexchange sites to restrict your post to a single question and do some preliminary research via google or your search engine of choice. But I don't see any harm in answering it here nonetheless.

The real answer is that naming conventions for genes/proteins are incredibly inconsistent and counterintuitive. You can't really draw any conclusions about a gene/protein based on its name alone. For example, some proteins such as p53 are named for biochemical features, others are named for phenotypes observed when their gene is mutated, such as for instance the gene decapentaplegic in fruit flies. Often they involve little jokes, such as the gene called mothers against decapentaplegic which represses the decapentaplegic mutant phenotype. So it's all a bit of a mess and people kind of assume that you'll follow along, which can be confusing.

It's a little rarer that gene names differ substantially from the names of the proteins they encode, but not unheard-of. This is indeed the case with your TP53/p53 example. See a short explainer here or a longer one here. What you'll notice is that there are in fact multiple different names for the same gene/protein. What you'll often see is that the shorter name (p53) is the one that people tend to use colloquially.

For your second question regarding how mutations are propagated, the short answer is that mutations are inherited from mother to daughter cells. Cells don't generally "infect" other cells with mutations (though in biology there are always exceptions). I think that what you are interested in is the difference between "somatic" and "germ-line" mutation.

Hope that helps.


How do I carry out an ancestry/admixture test on a single VCF file?

A modified implementation of Vivek's answer.

peddy is a Python package that samples an input .vcf at ~25000 sites and projects onto a principal component space built on 2504 thousand genome samples. The author has extensive documentation of the tool's features and a link to the preprint.

I downloaded the .vcf and .vcf.tbi for the NA12878 sample from Genome in a Bottle's ftp here. Then, created a custom .ped file, NA12878.ped, with the contents:

NA12878 HG001   0       0       2       0

At the command line:

$ peddy --plot --prefix myvcf HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz NA12878.ped

The output files all have the prefix myvcf., here's myvcf.pca_check.png


"Sequence Duplication Levels" module still fails after pre-processing Illumina data

FastQC assumes that all samples are for whole genome sequencing and will flag them as failed if they differ too much from that assumption. This will, for example, cause essentially all RNA-seq, ChIP-seq, and ATAC-seq samples to fail in one module or another. This is not any cause for concern and is completely expected. Primarily concern yourself with whether all of your samples are similar in their metrics.


I was using dudi_pca incorrectly. The supplied parameter to as.matrix() should have been dudi_pca$li:

rtsne_out <- Rtsne(as.matrix(dudi_pca$li), pca=FALSE)


Sophisticated algorithms in biological systems

There are many examples of biological systems doing computation, according to my understanding:

Regulatory networks encode decision trees of significant complexity: wikipedia.
Alan Turing did significant work on organismal morphogenesis showing how simple chemical patterns of inducers can lead to complex outcomes, "Turing patterns".
All of organismal development, cell division, etc. are tightly regulated procedures of gene regulation and morphogenesis.
The inspiration for neural networks is brains.


UCSC packages/tools ucsc-netchainsubset and ucsc-fatotwobit not working on GNU/Linux

Short answer: respect the case of the tool name

Longer story below

Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.

ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin
c_rehash         libpng-config      mysql_config    openssl  pngfix
libpng16-config  my_print_defaults  netChainSubset  perror   png-fix-itxt

Then I typed netChainSubset and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called netChainSubset but it seems it can be called with netchainsubsetor NETCHAINSUBset for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the bin of where the tool is installed.


How to get unique somatic mutations for each individual patients

If what you want is to split the main VCF file into 1 file per sample, you could use bcftools query and view commands.

A similar question was asked on biostars, adapting Jorge Amigo's solution to your situation, you might use:

for sample in $(bcftools query -l $combined_vcf)
do
    bcftools view -c1 -Ov -s $sample -o ${$sample.vcf} $file
done


Personally I think its late and we're all tired.

"FileNotFoundError: [Errno 2] No such file or directory: 'ls_orchid.fasta'"

The traceback is pointing to line 2, the SeqIO as you point out.
In my opinion, the script should be in the same directory as the fasta file. You know how to identify the path in bash (pwd or echo $PWD), although you must be using Windows (Anaconda has \ rather than /). Alternatively the file names has a typo. You know the
for seq_record in SeqIO.parse("~\path\ls_orchid.fasta", "fasta"): # where path is the dir(s) leading to ls_orchid.fasta, but obviously use / if its Linux.

I would assume you can alternatively dump the .ipynb file in the location where the fasta file is.
As personal choice I prefer a formal IDE such as 'Visual Studio Code' or "PyCharm". In my personal opinion, I am aware some will disagree, Jupyter notebook good is for running code on a remote commercial server for data science work. However, the traceback - which you clearly have activated - looks good. I can see that you're Anaconda3 installation looks good, you must be using Windows (which I've never worked with) and you're running Jupiter from Anaconda3 ... so cool.


From the page you cited:

  During first-strand synthesis, upon reaching the 5’ end of the RNA template, the terminal transferase activity of the MMLV reverse transcriptase adds a few additional nucleotide

That usually are CCC (1st is almost always C and the other two most of the times), see: Zajac et al., 2013

These CCCs are used for priming the Template Switching Oligo (TSO), that has GGG and whatever you want to add to the sequence (for instance UMIs and adapters).

So, answering to your question the polimerase does not switch strands, it switches from mRNA as a template to TSO as template. Because of that is called template switching.

You can see it in the graph from the page you cited. In the template switching step polimerase switch from the mRNA to the TSO and continue with the synthesis of the first cDNA strand using TSO as template.


Using multidimensional scaling to visualize protein sequences by functionality

Some of the first MSA analysis, eg of codon usage of bacteria were performed using matrix factorization (see old papers from Des Higgins). His group also published discriminative methods for finding functional residues using supervised correspondence analysis

  Wallace IM1, Higgins DG BMC Bioinformatics. 2007 Apr 23;8:135.
  Supervised multivariate analysis of sequence groups to identify
  specificity determining residues.

The method is available in the R package ade4 and Bioconductor package made4

Also have a look at the phylogenetic R package adephylo and several related packages


Separation of mixed plasmid DNA sequences post whole-plasmid sequencing

You would expect to have high coverage, given the plasmids are short, so de novo assembly would be likely very easy.  Given that each plasmid is present in different multiples, you would expect different coverage on each plasmid, so it might be best to approach it as a metagenome-type or transcriptome assembly, rather than a classic genome assembly.

Alternatively, Google search just told me there is a tool called plasmidSPADES designed for just your purpose.  SPADES itself is a reliable genome assembler, so I would imagine plasmidSPADES to be a good choice.  It does assume that there is also some whole bacterial genomes in the sequencing mix though.

As @terdon commented, it does depend on whether you have many identical copies of each plasmid, or many similar but not identical plasmids in the mix.


A contact map predictor will usually output the probability of a contact for every residue pair in a protein, often ommitting very close pairs (e.g. less than 6 residue separation).
To calculate the L/n score, rank the predicted contacts descending by probability and look at the first L/n contacts. The fraction of these that are actually contacts in the native structure is the L/n contact precision score.
Often, you would filter the list to look at only contacts of a certain sequence separation, e.g. long range contacts between residues of >23 separation (definitions vary between papers). If you have a good predictor this score is high for L/10 and falls as you go to L/1, which is expected as proteins don't usually have that many long range contacts to find.


snakemake - using output file not explicitly declared as input for subsequent rule

Another way is to add one command after minimap: touch done.txt. This creates an empty file after minimap is finished. This can then be used in output rule for snakemake.

But specifying all the output files may not be a bad idea if these outputs are used in subsequent rules.


How to transform a DNAStringSet from the Bioconductor package Biostrings to a data frame?

According to the documentation (?Biostrings::DNAStringSet):

  width(x): A vector of non-negative integers containing the number of
            letters for each element in x.  Note that width(x) is
            also defined for a character vector with no NAs and is
            equivalent to nchar(x, type="bytes").

  names(x): NULL or a character vector of the same length as x
            containing a short user-provided description or comment for
            each element in x.  These are the only data in an
            XStringSet object that can safely be changed by the user. All
            the other data are immutable!  As a general recommendation,
            the user should never try to modify an object by accessing
            its slots directly.

  as.character(x, use.names=TRUE): Converts x to a character vector
            of the same length as x.  The use.names argument controls
            whether or not names(x) should be propagated to the names
            of the returned vector.

This should be enough for you to create a data.frame:

dss2df <- function(dss) data.frame(width=width(dss), seq=as.character(dss), names=names(dss))

Also, if you have a huge DNAStringSet, perhaps you do not want to convert all the sequences to characters; storing sequences as characters supposedly is less efficient. For plotting it might be sufficient to have just widths.

qplot(width(dss), geom='histogram'))


How can I bypass 'out of memory' problem using a HISAT2 for human genome indexing?

You haven't allocated the requisite 160GB of RAM required. Or rather, you allocated 160GB of memory, which is what you were told to do, but they most likely meant 160GiB.

Computers use binary, not decimal, so a gigabyte is 1024 megabytes. This is the same annoying mix-up people run into when their hard drives don't have the space the marketing department says it does. By requesting 160000M in your job script, you actually requested a few gigabytes less than what your source says you need.

$$
160000 \text{ megabytes} \times \frac{1 \text{ gigabyte}}{1024 \text{ megabytes}} = 156.25 \text{ gigabytes}
$$

$$
160 \text{ gigabytes} \times \frac{1024 \text{ megabytes}}{1 \text{ gigabyte}} = 163840 \text{ megabytes}
$$

If "160GB" is what you need, then 163840M is the correct value for the job script.

I would honestly suggest cutting the Gordian knot and not dealing with this "gigabytes vs. gibibytes" ordeal by simply using the G suffix.1

#SBATCH --mem=160G


Converting Ensembl Gene IDs to Entrez Gene IDs through biomart

The ID you need is the NCBI gene ID, which is the same as the EntrezGene ID.


Confidence Interval with Wilcoxon Test in Python for log-normal Distribution


Which values from Seurat::FetchData function are to be used for correlation analysis between genes?

You should use normalized data. The function has a parameter slot="data" that let's you select to access the normalized data.
Any reason for rounding to 1 decimal place here? I imagine that would have a substantial impact on the correlation you compute.


GOseq analysis with evidence code filter

A bit weird to answer my own question, but for who's interested I can show the 'hard' way using biomaRt and then use goseq with gene2cat defined.

So first you need biomaRt to get annotation of the genes. I have ensembl genes from mm10.

library(biomaRt)

table <- read.table(file, header = T, stringsAsFactors = F)

ensembl <- useMart("ENSEMBL_MART_ENSEMBL")

ensembl <- useDataset("mmusculus_gene_ensembl", mart = ensembl)

allGenes <- table$Geneid

GOs <- getBM(attributes=c('ensembl_gene_id', 'go_id', 'name_1006', 'go_linkage_type', 'namespace_1003'),
      filters = 'ensembl_gene_id',
      values = allGenes,
      mart = ensembl)

dim(GOs)
[1] 477082      5

table(GOs$go_linkage_type)
          EXP    IBA     IC    IDA    IEA    IEP    IGI    IMP    IPI    ISA    ISM    ISO    ISS    NAS     ND    RCA    TAS
 36931     36  27918    525  34038 210114    630   6891  30219   6205   4234     20  77365  24673    609  12821     99   3754

# I only want BP and no IEA evidence code
GOs_no_IEA <- GOs[GOs$go_linkage_type != "IEA" & GOs$namespace_1003 == "biological_process",]

dim(GOs_no_IEA)
[1] 111944      5

Then you can continue with goseq.

library(goseq)

All <- read.table(input_All, sep = "\t", header = T)

All_genes <- as.character(All$Geneid)

TopTable <- read.table(input, sep = "\t", header = T)

sig_genes <- as.character(TopTable$Ensembl)

# Make gene vector

gene.vector <- as.integer(All_genes %in% sig_genes)

names(gene.vector) <- All_genes

gene.lengths <- All$Length

pwf <- nullp(gene.vector, "mm10", "ensGene", bias.data = gene.lengths)
head(pwf)

GO.wall <- goseq(pwf, "mm10", "ensGene", gene2cat = GOs_no_IEA[, c(1:2)], use_genes_without_cat = F, test.cats = "GO:BP")


How to calculate mutation rate and mutation sites in a genome using FASTA file?

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


How to do BLASTP for short sequences (<20 aa) effectively?

You should increase the -evalue 10 default option to a higher value because this value it's adjusted to the sequence database size. Check this table, maybe this can help you to decide which parameters are the best option for you.

    OPTION      TASK          DEFAUL    DESCRIPTION
    word_size   blastp        3         Word size of initial match.
    word_size   blastp-short  2         Word size of initial match.
    gapopen     blastp        1         Cost to open a gap.
    gapopen     blastp-short  9         Cost to open a gap.
    gapextend   blastp        1         Cost to extend a gap.
    gapextend   blastp-short  1         Cost to extend a gap.
    matrix      blastp        BLOSUM62  Scoring matrix name.
    matrix      blastp-short  PAM30     Scoring matrix name.
    threshold   blastp        11        Minimum score to add a word to the BLAST lookup table.
    threshold   blastp-short  16        Minimum score to add a word to the BLAST lookup table.
    window_size blastp        40        Multiple hits window size, use 0 to specify 1-hit algorithm.
    window_size blastp-short  15        Multiple hits window size, use 0 to specify 1-hit algorithm.


Aligned base strand in pysam `pileupcolumn`

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


Why most aligners do not output the "X" CIGAR operation?

The SAM format originally had only M, I, D, N, S, H, and P CIGAR operators.  See the original SAM specification (if you can view Apple Pages documents) and Table 1 in The Sequence Alignment/Map format and SAMtools (Li et al, 2009).  This was in line with previous tools using CIGAR strings, notably exonerate which introduced them with just the M, I, and D operators.

BWA-backtrack was written contemporaneously with the SAM format in 2008 and 2009 (and published in May 2009).  Its ChangeLog shows that it was outputting mismatch information in an MD tag from January 2009 and that MD was defined in the SAM specification of the time (and that the tag value's syntax was somewhat in flux in February 2009).  The MD tag is also described in that early v1.0 Pages-formatted SAM specification.

The = and X CIGAR operators were introduced later in SAM v1.3 as a result of this lengthy samtools-devel mailing list thread.  The characters used for the operators and the initial implementations were essentially in place by November 2009.

Since then, the =/X operators have not really taken over from M, as you've seen.  There's probably a number of factors contributing to this:

Their later introduction, long after the basic M/I/D operators and MD tag were well-established;
Conceivably their being specific to SAM and unavailable in other CIGAR flavours;
The MD tag still provides more information — X doesn't tell you what the mismatched reference bases were.


Why is this makeblastdb command not working?

This looks like a bug in makeblastdb. Removing the final | | | from your sequence's description makes it work:

>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15

So does removing everything after the 1st space:

>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48

There is nothing in the FASTA definition that would explain this and, even more weirdly, the presence of the | | | alone isn't enough to trigger the issue. I played around with this a little and it looks like it doesn't like the combination of + or -along with empty | |, but there's no clear pattern I can see:

These fail:

>| | |15-| | |
>| | |15+| | |
>| | |15-1| | |
>| | |15+1| | |

These work:

>| |15+| |
>| | |151| | |

So yeah, sounds like a bug. I suggest you report it to the NCBI devs. In the meantime, a simple fix would be to just keep the first space-separated field of your fasta header. Pass your fasta file through this:

awk '{print $1}' fasta.fa > new.fa

That will print the 1st field of each line and, assuming you don't have spaces in your sequence, it should effectively sort this out for you.


Pseudo-temporal ordering in heterogeneous populations

Regarding the first part of the question: scRNA-seq is a rapidly developing field so it may be hard to talk about "widely accepted tool for doing pseudo-temporal ordering from scRNAseq data". Few of the tools aiming to do just that include Monocole, Waterfall or Sincell (see this paper for references https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4122333/)

The second part is a bit more complex. Many aspects have come together here, e.g. cell-cycles, sub-populaitons identification and pseudo-temporal ordering, to get a true reflection of the biological processes. There are efforts on all (and more) these fronts (see again e.g. the above paper, though not the latest) and there are most likely people working on their integration. I'm not aware of any studies published yet at this depth though


Oxford Nanopore mapping Quality and sequencing error

Where i know with the phred score mapQ of 60 means 99.9999% of the bases are called correctly.

That's not at all what the MAPQ means, it means that the sequence is very likely to have originated from that place in the genome, whether the underlying sequence is correct is another question all together.

For long nanopore reads, sequence quality is not well correlated with MAPQ. MAPQ will tend to be either very low (0 or 1) or very high (60) since you tend to have kb+ reads and those will largely be uniquely assignable in the genome. The sequencing quality of nanopore reads is only an issue when you care about something at the base level, such as SNP calling or finding exact splice junctions. Even then, you mostly need sufficient coverage to get around this.


select unique names while removing the duplicates from a column

That sound like a job for a perl one-liner:

perl -ane 'unless(exists $seen{$F[0]}){print $_; $seen{$F[0]}++}'

Explanation:
The input is read line-wise and printed back to the command line if the ID in the first column has not been seen before.

perl -ane ' reads input line-wise and splits fields to array @F
unless(exists $seen{$F[0]}){ checks if ID exists in hash %seen and continues if that is not the case
print $_; prints the original input to STDOUT
$seen{$F[0]}++}' adds the ID to the hash %seen


How to find GEO data sets using Drug Bank ID and bioDBnet

You could take a look at the rentrez package for R. You can use the entrez_search function to get the GEO IDs for the studies related to your query, and then the entrez_summary function for the results. So, for example:

x <- rentrez::entrez_search(db="gds", term=paste0("tolvaptan AND gse[ETYP]"),
                            use_history = T)
res <- rentrez::entrez_summary(db="gds", id=x$ids)
str(res)

Where res will be a list containing objects of the esummary class. Then, you can use the following line to get the GSE accession numbers related to tolvaptan:

lapply(res, `[[`, "gse")

So you could put this function in a loop replacing "tolvaptan" with your drugs:

drugs <- c("drugA", "drugB", "drugC")
for(drug in drugs){
    x <- rentrez::entrez_search(db="gds", term=paste0(drug, " AND gse[ETYP]"))
    res <- rentrez::entrez_summary(db="gds", id=x$ids)
    lapply(res, `[[`, "gse")
}

Take a look at the rentrez documentation for more details.


No output when trying to obtain protein sequence from PDB file

That PDB entry is an RNA molecule, so there is no protein sequence to extract.

If it extracted the RNA sequence by default here you can see there could be an ambiguity as to whether it was an RNA sequence or a protein sequence, as the alphabets share characters.


How to get Uniprot and swissprot cross references using Biopython

I believe the attribute you should be looking for is dbxrefs:

[x for x in record.dbxrefs if x.startswith('GO:')]

Apparently SeqIO handles uniprot text files differently than the SwissProt module:

from Bio import SwissProt
record = SwissProt.read(open('single_record_from_uniprot.txt'))
record.cross_references

So you should change you code to use the SwissProt module to get the behavior you expect


What is the actual cause of excessive zeroes in single cell RNA-seq data? Is it PCR?

It may be necessary to distinguish between methods that use unique molecular identifiers (UMIs), such as 10X's Chromium, Drop-seq, etc, and non-UMI methods, such as SMRT-seq.  At least for UMI-based methods, the alternative perspective, that there is no significant zero-inflation in scRNA-seq, is also advocated in the single-cell research community.  The argument is straight-forward: the empirical mean expression vs. dropout rate curve matches the theoretically predicted one, given the current levels of capture efficiency.

Examples

Svensson Blog

A couple of blog posts from Valentine Svensson argue this point rather pedagogically, and include citations from across the literature:

Droplet scRNA-seq is not zero-inflated

Count-depth variation makes Poisson scRNA-seq data negative binomial

bayNorm

There is a more extensive preprint by Tang, Shahrezaei, et al. (BioRxiv, 2018) that claims to show a binomial model is sufficient to account for the observed dropout noise.  Here is a snippet of a relevant conclusion:

  Importantly, as bayNorm recovered dropout rates successfully in both UMI-based and non-UMI protocols without the need for specific assumptions, we conclude that invoking zero-inflation models is not required to describe scRNA-seq data.  Consistent with this, the differences in mean expression levels of lowly expressed genes observed between bulk and scRNA-seq data, which were suggested to be indicative of zero-inflation, were recovered by our simulated data using the binomial model only.

Multinomial Modeling

There is also a very clearly written preprint by Townes, Irizarry, et al. (BioRxiv, 2019) where the authors consider scRNA-seq as a proper compositional sampling (i.e., multinomial process) and they come to a similar conclusion, though specifically for UMI-based methods. From the paper:

  The multinomial model makes two predictions which we verified using negative control data. First, the fraction of zeros in a sample (cell or droplet) is inversely related to the total number of UMIs in that sample. Second, the probability of an endogenous gene or ERCC spike-in having zero counts is a decreasing function of its mean expression (equations provided in Methods). Both of these predictions were validated by the negative control data (Figure 1). In particular, the empirical probability of a gene being zero across droplets was well calibrated to the theoretical prediction based on the multinomial model.  This also demonstrates that UMI counts are not zero inflated.

Furthermore, by comparing raw read counts (prior to UMI-based deduplication) and UMI counts, they conclude that PCR is indeed the cause of zero-inflation:

  The results suggest that while read counts appear zero-inflated and multimodal, UMI counts follow a discrete distribution with no zero inflation (Figure S1). The apparent zero inflation in read counts is a result of PCR duplicates.

I highly recommend giving this a read, especially because it nicely situates other common generative models (e.g., binomial, Poisson) as valid simplifying assumptions of the multinomial model.

It should be noted that this same group previously published a work (Hicks, Irizarry, et al. 2018), mostly focused on non-UMI-based datasets (SMRT-seq), where they showed evidence that, relative to bulk RNA-seq, there was significant zero-inflation.


Compile with gcc or remove inline from librna/fold.c:776 and librna/fold.c:846.


Restricting match output for multiple sequence alignment using mafft?

Q - Is there some parameter for this?

No, there is no parameter for this. Getting the conserved regions is something that you will have to perform post-alignment.

How to get the conserved regions?

You will need to define what exactly you mean by conservation and calculate a conservation score for each position of the alignment.

For instance you can look at each column of the alignment and calculate the percent of queries that had the same base/AA at that position. If using AAs you can define conservation more loosely depending on what properties are important to you, or using a scoring matrix.

Once you have a co-occurrence matrix you can define the regions of interest based on a conservation threshold that you define.

Is there software to find conserved columns?

There is a program MstatX which outputs conservation scores of individual columns.

How MstatX works

The example of MstatX provides this example alignment. Each column is progressively less conserved than the last. The last column is conserved, but has some gaps -- this lowers the conservation score.

   col          1
 seq1  1234567890
 seq2  DDDDDDIPDL
 seq3  DDDDDDIPVL
 seq4  DDDDDDIPYL
 seq5  DDDDDDIPA-
 seq6  DDDDDDLWT-
 seq7  DDEDEELWK-
 seq8  DDEDEELWP-
 seq9  DDEDEFVSR-
 seq10 DEEFFFVSH-

The conservation scores for each column were:

1       2.71272e-08
2       0.243433
3       0.312789
4       0.243433
5       0.483364
6       0.480391
7       0.498757
8       0.498757
9       0.973514
10      0.267319

How to get conserved regions

The python program below will output all windows of conserved regions when using a fasta alignment and the default output from MstatX.

Positional args

1. MstatX_alignment_example.mali, is the fasta alignment file.
2. output.txt is the default output of mstatx when called with mstatx -i MstatX_alignment_example.mali.
3. 0.4 is the minimum conservation score
4. 3 is the minimum number of bases that must have that conservation score.

Calling this program with the following command on the above example:

python print_conserved_regions.py MstatX_alignment_example.mali output.txt 0.4 3 prints:

{1,4}

Meaning that there was only one window at least 3 bases long in which all positions had a conservation score below 0.4. This window was positions 1 through 4 in the alignment.

#call this file cons_windows.py
import sys
from Bio import AlignIO

align_filename = sys.argv[1]
cons_filename = sys.argv[2]
cutoff = float(sys.argv[3])
minchar = int(sys.argv[4])

def parse_conservation(filename):
    cons = {}
    with open(filename, "r") as f:
        for line in f:
            if line.strip():
                temp = line.split()
                cons[int(temp[0])] = float(temp[1])
    return cons

align = AlignIO.read(align_filename, "fasta")
cons = parse_conservation(sys.argv[2])

i = 1
while i < max(cons):
    if cons[i] <= cutoff:
        j = i
        done = False
        while j < max(cons) and not done:
            if (cons[j] > cutoff) and (j-i >= minchar):
                print("{{{},{}}}".format(i,j-1))
                done = True
                i = j
            j+=1
    i+=1


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

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

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


Consensus sequence from SAM or BAM file?

A few years ago I used AlignGraph (paper and code) for a task like this. I gave AlignGraph my raw reads, so I'm not sure if it accepts SAM/BAM input or if you have to go back to the original Fastq reads.

I'm not sure if any related/improved tools have been published since AlignGraph. Internet searches for "reference-guided assembly" return mostly transcriptome assembly tools. It might be worth looking for articles that cite the AlignGraph paper.


I think these definitions still have a degree of subjectivity and I don't think they are mutually exclusive.

Convergent amplification basically means that "similar" cancers/neoplasms have "similar" amplification patterns, i.e. at around the same loci (see Myllykangas S et al.,2006 for an example).

Focal genome amplification is defined as "CNA [copy number alteration] of limited size, frequently enriched for cancer driver genes" (Krijgsmanthis et al., 2014).


Variants in Exome that affect miRNAs related to Schizophrenia

I would run them through the Ensembl VEP, which will tell you which genes they hit and the effect they have on them. Then I would filter by consequence "mature_miRNA_variant", which will find you all the variants in your file that fall in the mature miRNAs. Or you could filter instead by the gene biotype "miRNA" which will tell you if the variants fall anywhere in the miRNAs (which will include the full stem loop structure). I would then export the list of miRNAs and try to find out about their targets using miRbase.


How to interpret the SMRT Link base modification algorithm output?

It sounds like you're looking for the modifications.csv file, rather than the modifications.gff file. The .gff file is for viewing in SMRT View (or other viewers). You can see schema for both files at PacBio's Methylome Analysis Technical Note. Specifically, modified_base means the algorithm could not determine the exact base and modification.

I also found a video tutorial for performing the analysis in SMRTLink on YouTube

EDIT: Per PacBio's Methylome Analysis Technical Note documentation, the .csv "...is also produced when motif analysis is not active, such as with the RS_Modification_Detection analysis protocol." Not sure why it's missing, perhaps because the algorithm couldn't clearly identify any of the modifications. You'll probably need to contact tech support on this issue.

In either case, the .gff you posted includes the significant hits (p < 0.01), but the algorithm was unable to classify them, likely because of insufficient coverage. PacBio's Methylome Analysis Technical Note documentation gives specific coverage recommendations for each event type, but 100x is the lower bound. Lots of good information there.

Probably not what you wanted to hear, but I hope that gets you going.


Bruker MALDI-TOF bacteria species identification scoring algorithm

For identification scoring, the software correlates signal
intensities of matched signals of mass spectra. The three scores
obtained from such a procedure are multiplied and normalized to a
value of 1,000 and its common logarithm. Log scores over 2 are
considered as reliable identification of a bacterial species (Table
4), whereas log scores over 1.7 generally indicate reliable
identification of bacterial genera. Log scores of 3 are obtained when
spectra are matched with themselves.

Sauer, S., Freiwald, A., Maier, T., Kube, M., Reinhardt, R., Kostrzewa, M., & Geider, K. (2008). Classification and identification of bacteria by mass spectrometry and
computational analysis. PloS one, 3(7), e2843.
Going further, from a random thread on ResearchGate:
https://www.researchgate.net/post/Does_anyone_know_how_is_the_p-value_calculated_by_Bruker_Biotyper_during_MALDI-MS_fingerprinting_of_bacteria

Number of peaks matched in your unknown sample divided by the number
of peaks in the unknown sample, this is given a score (out of 10)
Number of peaks matched in your unknown sample divided by the number
of peaks in the reference sample, this is given a score (out of 10)
Score of relative intensities of matching peaks (out of 10)
Total score out of 1000
Log10 score out of 3
There is no relationship of the score with a p-value nor you can
calculate the statistics because the bruker database is not open and
you can not know the values for what the score was calculated to do
the statistical analysis


Get the mapping statistics of a single read from a BAM file

samtools stats seems to be able to do most of this, excluding the CIGAR-string parsing stuff (i.e. INDELs):

$ samtools view -h mapped.bam | grep -e '^@' -e 'readName' |
    samtools stats | grep '^SN' | cut -f 2-

raw total sequences:    2
filtered sequences: 0
sequences:  2
is sorted:  1
1st fragments:  2
last fragments: 0
reads mapped:   2
reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
reads unmapped: 0
reads properly paired:  0   # proper-pair bit set
reads paired:   0   # paired-end technology bit set
reads duplicated:   0   # PCR or optical duplicate bit set
reads MQ0:  0   # mapped and MQ=0
reads QC failed:    0
non-primary alignments: 0
total length:   13632   # ignores clipping
bases mapped:   13632   # ignores clipping
bases mapped (cigar):   13502   # more accurate
bases trimmed:  0
bases duplicated:   0
mismatches: 3670    # from NM fields
error rate: 2.718116e-01    # mismatches / bases mapped (cigar)
average length: 6816
maximum length: 6816
average quality:    8.9
insert size average:    0.0
insert size standard deviation: 0.0
inward oriented pairs:  0
outward oriented pairs: 0
pairs with other orientation:   0
pairs on different chromosomes: 0

For a length filter, an awk length check can be slotted in instead of the grep:

$ samtools view -h mapped.bam |
    awk -F'\t' '{if((/^@/) || (length($10)>500)){print $0}}' |
    samtools stats | grep '^SN' | cut -f 2-

raw total sequences:    131
filtered sequences: 0
sequences:  131
is sorted:  1
1st fragments:  131
last fragments: 0
reads mapped:   131
reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
reads unmapped: 0
reads properly paired:  0   # proper-pair bit set
reads paired:   0   # paired-end technology bit set
reads duplicated:   0   # PCR or optical duplicate bit set
reads MQ0:  0   # mapped and MQ=0
reads QC failed:    0
non-primary alignments: 0
total length:   569425  # ignores clipping
bases mapped:   569425  # ignores clipping
bases mapped (cigar):   453884  # more accurate
bases trimmed:  0
bases duplicated:   0
mismatches: 113755  # from NM fields
error rate: 2.506257e-01    # mismatches / bases mapped (cigar)
average length: 4346
maximum length: 16909
average quality:    9.1
insert size average:    0.0
insert size standard deviation: 0.0
inward oriented pairs:  0
outward oriented pairs: 0
pairs with other orientation:   0
pairs on different chromosomes: 0


How to count the number of mapped read in 100-bp window from a BAM/SAM file

one-liner

Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ($4 >=1) and the 500 in ($4 <=500) to set your window. Change "hg19" to your target sequence.

Note, this one-liner does not double-count reads because of uniq.

samtools view file_sorted.bam | \
    awk '{ if (($3 == "hg19") && ($4 >=1) && ($4 <=500)) { \
  print($1) }}' | sort | uniq | wc -l

python 3

If you need something a little more programmatic to get a known region you should consider trying the pysam package.

import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
region_set = set()
start = 100
stop =  500
counter = 0
for read in samfile.fetch('interesting_contig', 100, 500):
    region_set.add(read.query_name)
    counter += 1
print("{} unique reads in [{}, {}]".format(len(region_set), start, stop))
print("{} alignments in [{}, {}]".format(counter, start, stop))

For my bam file, this program outputs:

1618 unique reads in [100, 500]
2487 alignments in [100, 500]

python 3, non-overlapping windows

It is important to note that the samfile.fetch() function conforms to 1-based indexing, consistent with SAM files, unlike pysam's otherwise default 0-based indexing that matches the BAM-indexing.

import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
print(samfile.references)
print(samfile.lengths)
window = 500
print("ref\tstart\tstop\tuniq_reads")
for i in range(len(samfile.references)):
    refname = samfile.references[i]
    seqlen = samfile.lengths[i]
    for j in range(1, seqlen, window):
        stop = j+window-1 if j+window-1 < samfile.lengths[i] else samfile.lengths[i]
        region_set = set()
        for read in samfile.fetch(refname, j, stop):
            region_set.add(read.query_name)
        print("{}\t{}\t{}\t{}".format(refname, j, stop, len(region_set)))

For the small bam file I looked at the output was:

ref         start  stop  uniq_reads
OdonB       0      500   1668
OdonB       501    1000  2647
OdonB       1001   1500  129400
OdonB       1501   2000  598
OdonB       2001   2500  538
OdonB       2501   3000  872
OdonB       3001   3500  561
OdonB       3501   4000  49120
OdonB       4001   4500  48433
OdonB       4501   4806  28121


using vcffilterjdk:

you want to keep the variant if there is no INFO/CAF.

Otherwise, get All values of CAF as a String list (because it is declared as a string in the VCF header...), stream those data, skip the first value for REF , convert each string to double, you want to keep the variant if any allele is rare (< 0.01)

java -jar dist/vcffilterjdk.jar -e 'return !variant.hasAttribute("CAF") || variant.variant.getAttributeAsStringList("CAF","1.0").stream().skip(1L).map(A->Double.parseDouble(A)).anyMatch(A->A<0.01);' input.vcf


Why doesn't an (Entrez eutils) einfo request for "gene" return the link gene_nucleotide or gene_nucleotide_pos links?

I've been corresponding with NLM about this issue and I finally took the time to try out their suggestion (which personally I found hard to see between the lines and which is not a discrete solution, but rather a very time-consuming manual process containing false positives because they say that to get a formal and comprehensive response to my query:

Collecting the linkname will be a difficult task, that will take time and coordinate/check with relevant groups maintaining individual databases. Your patience will be greatly apprecired. [sic]

It took me a bit to decipher their suggestion on how to manually find all links, so I will share what I've learned.  There is (currently) no codified or formal means to obtain all possible link names and I infer that many links exist simply because they are utilized in the internal workings of their website.  I could be wrong about that, but suffice it to say that there are many undocumented links and they change constantly.
You can get a list of filter items using each individual database's advanced search web interface, "most of which should be the linkname from that source database to other target database".  (So this should result in all possible links among a series of false positives.)
Here's how you do that.  Let's take the gene database as an example:

Go to the advanced search web page for that database, e.g. https://www.ncbi.nlm.nih.gov/gene/advanced
In the "Builder" section, select "Filter"
Click "Show Index List" to the right of the search term field next to the select list where you selected "Filter" (you will see a multi-select list appear below the search field)
Click "Next 200" until you find a selection in that new select list matching the database name followed by "all" and a number in parentheses, e.g. gene all (29994947) and click to select that item
Click "Refresh index" at the lower right corner next to the multi-select list

The multi-select list will be repopulated and most of the items will be linknames (though you have to replace spaces with underscores in order for those linknames to work with the elink utility).
Doing this for the gene database and then scrolling through just links to nuccore (and its alias "nucleotide" - ignoring est and gss), you will find:

gene nuccore
gene nuccore clust
gene nuccore mgc
gene nuccore pos
gene nuccore refseqgene
gene nuccore refseqrna
gene nucleotide
gene nucleotide clust
gene nucleotide mgc
gene nucleotide mgc url
gene nucleotide pos

If you look only at the linknames in the documentation or via einfo (by supplying the gene database and looking at similar links to nuccore), you only get these links:

gene_nuccore
gene_nuccore_mgc
gene_nuccore_refseqgene
gene_nuccore_refseqrna
gene_nucleotide_mgc_url

So I believe that given this information and the ever-changing nature of these links, I believe I will allow the user the option to enter a linkname manually, if the set in the select list I generate via einfo does not contain the link they need.


Splitting fasta file into smaller files based on header pattern


RepeatMasker comes with a buildSummary.pl script.

You can input the .out file from your RepeatMasker run and there is a flag to enter your species

e.g.

buildSummary.pl -species Vertebrates rmrun.out > mysummary.tbl

In the mysummary.tbl, the summary will not assume your run to be Homo (default).


Proof of Breakpoint Reversal Sorting Approximation Algorithm

Denote the strips in your permutation as $\langle S_1, S_2, \cdots, S_n \rangle$. Suppose $k-1$ is in segment $S_i$, and $k$ is in segment $S_j$, for some $i, j \in [n], i \neq j$.

Recall by definition of $k$ that $k-1$ must be the maximum and therefore rightmost element in an increasing strip, and $k$ the minimum and therefore also rightmost in a decreasing strip. Consider the two cases:

$i < j$. The increasing strip comes before the decreasing. $k-1$ is on the 'right' end of $S_i$, but $k$ is on the 'right' end of $S_j$ also. In order for $k-1$ and $k$ to be adjacent, reverse the segments $S_{i+1}, \cdots, S_j$.
$j < i$. The decreasing strip comes before the increasing. $k$ is on the 'right' end of $S_j$, but $k-1$ is on the 'right' end of $S_i$ also. Again, we reverse $S_{i+1}, \cdots, S_j$.

So, "reversing the segment between $k$ and $k-1$" means to reverse all the strips in-between the two strips containing $k$ and $k-1$, excluding the 'leftmost' strip, but including the 'rightmost' strip.


How to annotate gene length to a list of gene symbols using UCSC data?

You could do a few command-line operations to answer this question. This assumes the use of hg38 assembly.
First, get a list of genes from GENCODE:
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

Then use grep to filter those genes with your list of HGNC symbols:
$ grep -wFf hgnc_symbols.txt genes.bed > filtered_genes.bed

You can run this through awk to get lengths:
$ awk -vFS="\t" '{ print $3-$2 }' filtered_genes.bed > filtered_gene_lengths.txt


Using R to import specific records from a PDB file

If your file is in PDB format as the one shown here in .pdb format then you can simply extract your desired fields using sed which is a programming language used to parse and transform text. As for this structure i.e 6C4A in .pdb format HETNAM and HETSYN can be extracted using sed as follows:

sed  -n -e '/HETNAM/ p' -e '/HETSYN/ p'  6C4A.pdb

HETNAM     EJA S-[(1Z)-2-CARBOXY-N-HYDROXYETHANIMIDOYL]-L-CYSTEINE
HETNAM     PYR PYRUVIC ACID
HETNAM      MG MAGNESIUM ION
HETNAM     3NP 3-NITROPROPANOIC ACID
HETNAM     PEG DI(HYDROXYETHYL)ETHER
HETNAM     EPE 4-(2-HYDROXYETHYL)-1-PIPERAZINE ETHANESULFONIC ACID
HETSYN     EPE HEPES


Are there any databases of templates for common bioinformatic file formats?

You mention Biopython, which contains tests: https://github.com/biopython/biopython/tree/master/Tests.

Some of the tests consist in reading files present in the folders listed in the above link. These files could be a starting point for a database of test files. Whenever one comes across a test case not covered with these files, one could construct a new test file and contribute it to Biopython, along with a test, or at least file an issue: https://github.com/biopython/biopython/issues.

That would be a way to contribute to Biopython while constituting a database of test files.


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

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

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


Resolution parameter in Seurat's FindClusters function for larger cell numbers

Assuming you have an informative selection of variable genes from which you have constructed a number of useful PCs, I'd run a number of iterations with FindClusters() as described in the other answer, then choose a level which overclusters the dataset (for example, clusters that are visibly separate on a t-SNE or other dimensionality reduction plot should definitely have their own number):

seuratobject <- SetAllIdent(seuratobject, id='chosen.resolution')

Then run:

Seurat::BuildClusterTree()
Seurat::FindAllMarkersNode()

Assessing the cluster markers for each node will hopefully give you a good idea on which clusters should be combined. Then you can "combine" the clusters and re-label the cells using something like:

library(plyr)
cell.labels <- seuratobject@ident
cell.labels <- mapvalues(cell.labels,
  from=0:16, # cluster numbers
  to=c('A', 'B', 'C', 'C', 'D', 'E', 'E', ... )) # etc

seuratobject <- AddMetaData(seuratobject, cell.labels, 'Combined.clusters')

The usefulness of the clustering will very much depend on the selection of variable genes, therefore, depending on the (diversity of the) dataset, you will want to experiment with selection parameters or subset the dataset and repeat the above procedure.


What is the perferred method of optimization or energy minimization of small molecules downloaded from PubChem?

The difference between different force-fields is not going to be major, it is the side steps which are.

When

If you are starting from a SMILES string, optimisation is a must obviously.

If you are using a 3D conformer from PubChem or even an actual sub-1 Å crystal structure from CSD, optimisation is nice for consistency.

Which

MMF94 is a solid choice. RDkit offers MMF94 (Merck force field 1994) and UFF (universal force field). The latter has some extra features, such as not kicking the bucket at atoms that aren't your usual suspects. ChemDraw offers MM2 and MMF94 IIRC, the former is very old. g16 in Gaussian (paid software) uses Gaussian Electrostatic Model, which is polarisable and will solve partial charges in a nice way (which I'll get to in a minute). OpenBabel uses MMF94 under the hood. You can even optimise with macromolecular forcefields, such as AMBER and CHARMM.

Given that the AutoDock uses a flexible ligand it should not matter too much anyway.

What else?

Autodock does not require conformers by virtue of the flexible ligand business (the sampler does rototorsions). Other programs might if they use rigid ligands (e.g. Fred) or they use a mixed approach and use a conformer library for speed (e.g. Rosetta ligand_dock). Conformer diversity is actually a bigger issue as conformers that are only slightly contorted, often even less than kB (0.6 kcal/mol), may not be sampled.

The thing that messes things is incorrect protonation. OpenBabel can protonate at a given pH, RDKit cannot. Not an issue for ligands with drug screens (cf. Lipinski rule), but is for natural substrates. However, for the protein it is essential to make sure acidic and basic residue are correctly protonated and that histidine have the correct tautomer —there are tools, c.f. MOE, that can do it but by far the best approach is reading up on the mechanism.

The next thing is dodgy Gasteiger Marsili partial charges —for Autodock 4 (and most other full/hybrid physics-based docking programs), partial charges are important: after all PDBQT has a Q in it for charge! For prepare_ligand4.py you give it a mol2 file. A major consideration is what format you are using for your molecules. This sounds a very silly thing to point out, but you'd be surprised how many trip up. An mol file (or an sdf file, which is the same but multi-entry) does not have an official partial charge column (although they are often added as comments after M END). So whatever force-field you chose make sure that the program you use gave you partial charges are not formal charges. This does not quite apply to AutoDock Vina, which is charge independent and more reliant on machine learning derived factors.

Verdict

Open babel is nice in that it can start with a SMILES, handles multiple formats including mol2 and PDBQT, can protonate at pH 7, uses MMF94 and does Gasteiger Marsili charges (and its confab algorithm is very good, not need here though).


IonTorrent variant caller failing, lack of ZM: tags in BAM files

It turns out, after speaking to tech support, that these are old files from a previous version of the software so won't run through variant caller.

The files themselves need reanalysing from the initial FASTQ -> BAM with the newest version.


Is there any publicly available multi-omics dataset?

There is a database called OmicsDI, where one can search for multi-omics datasets.

Here's a link of the associated publication (Perez-Riverol, Yasset, et al. "Omics Discovery Index-Discovering and Linking Public Omics Datasets." bioRxiv (2016): 049205.) for more details.


scRNA-seq,10x cellranger pipelines,low custom tdtomato gene content! looking for help！

A few suggestions:

Manually inspecting the BAM file with a genome browser like IGV, as @GWW suggested and checking the qualities of the alignments in your BAM/SAM file (SAM is the human readable version of BAM, there are tools to convert one to another);
Using another single cell analysis tool like Seurat so that you can manually set thresholds for QC and not lose to many (tdtomato) expressing cells at the QC step (I don't know about your experiment but if transfected, cells could be at a worse condition than non-transfected ones and hence might not pass QC);
Double checking your tdtomato reference sequence, I don't know about tdtomato but for example GFP has flavors like EGFP, similar but different sequences, might explain low alignment score if you observe in step 1;
Taking the limitations of scRNA-seq (especially droplet-based) technologies into account. The figure below is from the CITE-seq paper (from supp info), Simultaneous epitope and transcriptome measurement in single cells (Stoeckius et al., 2017), showing that the protein - RNA measurements of a given gene might show poor correlation.


Create galaxy local toolshed using Planemo

--shed_target local will only work if you're running a Galaxy toolshed on your local computer, since it's reading the .shed.yaml file you made with planemo shed_init and trying to create a repository with that name or a local (in this case) server. Planemo itself won't start a Galaxy toolshed, since it's fairly unusual for people to want to do testing involving creating local repositories. What you probably want to do instead is use planemo test instead. Once the planemo testing passes then use the various shed* commands to upload your wrapper to the toolshed or testtoolshed.


How to align genomic sequence with corresponding amino acid sequence

Have a look at exonerate's protein2genome!
From the documentation:

  Aligning a protein to genomic sequence:

  Similarly, it is possible to align a protein sequence to the genome, (similar to GeneWise, but with
  heuristics).

  exonerate --model protein2genome query.fasta target.fasta

  This model will allow introns in the alignment, but also allow frameshifts, and
  exon phase changes when a codon is split by an intron.


e-utils api: get chromosome position for each snp on a gene

If your dataset is limited to human SNPs, and you know the list of SNPs, I would advise taking a look at the dedicated NCBI Variation Services.  This would be one request for each SNP, but if you fire off 1-2 per second, I think you'll be ok.
For example, for rs1490381980, you could try something like:
curl -X GET 'https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/1490381980' \
   | jq '.primary_snapshot_data.placements_with_allele[] | "\(.seq_id) \(.alleles[0].allele.spdi.position)"'

It would be better to query for the SNP and parse the JSON in python, but jq helps provide an idea of what is possible.
Alternatively, you could download the JSON or VCF data from the NCBI dbSNP FTP site, and process out the information you need.


remaining human genome variation that hasn't been sequenced?

In theory, almost any base in the human genome may mutate, so you have billions of variants to go. Ok, this is not so useful.

A related and potentially useful question is: given a human, what is the fraction of his/her variants seen in the the previously sequenced samples? For a Wright-Fisher population, there is an analytical answer: x% of variants on a haplotype has frequency below x%. Replace "x" with "1": 1% of variants on a haplotype have frequency below 1%. The 1000 genomes project concluded that 99% of SNPs at frequency 1% or higher have been called by the project. If we all came from a Wright-Fisher population, 99% of of SNPs have been called in 1000g. Given that the project also calls SNPs at lower frequency, nearly all SNPs in a newly sequenced sample are present in 1000g.

The above is a theoretical analysis. In practice, there are two major complications. First, we are not a Wright-Fisher population. Due to recent population expansion, each haplotype harbors more new mutations than the Wright-Fisher model would predict. Given a known population history, it is actually possible to numerically compute $f(x)$, the fraction of variants with frequency below x (for Wright-Fisher, $f(x)=x$). I don't have this result, though. Second, 1000g and many other genome projects didn't call SNPs in repetitive regions and missed a significant portion of high-frequency indels and structural variations. With Illumina, you can only call ~30% of long deletions callable with PacBio data, but we only have a dozen of PacBio human genomes in public. There are still a lot we don't know.


Ribosomal RNA QC quantification using SILVA

SILVA is intended for metagenomic samples and thus does not contain the human ribosomal subunits.  While the 5.8S subunit is in the ensembl, ncbi, and ucsc annotations, and there are numerous rRNA "biotype" (or "biomol") entries, the 5S, 18S, and 28S subunit sequences are absent.  Even downloading rRNA entries in the rmsk table from ucsc will not return all 4 major human subunits.

You can however (as I discovered) download all rDNA gene sequences (including those 4 subunits) from the nucleotide db on ncbi using the search term "biomol_rrna[PROP]".  This returns 240,240 coding sequences (or 30,775 gene sequences).  You can then use this as a reference for mapping, blasting, or as input to rRNA QC tools.

The following e-Utilities commands from ncbi will download the gene sequences:

esearch -db nuccore -query "biomol_rRNA[prop]" | efetch -format fasta > all_rrnas.fasta

This is not a common analysis.  Usually, rRNA reads can just be ignored and analyses are unaffected, which is likely why there is a paucity of records for these sequences in annotation files, but mapping to just these sequences could make for a simple QC measure as a gauge for rRNA contamination.  While the sequences may be present in the genome, they're not in the annotations and may even be masked as repetitive, since there are so many copies.  So you have to be careful to count unique reads and not mappings.


Why indexing sequence file with minimap2?

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

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


By reading this thread on seqanswers and by comparing the data to TCGA, I figured out

raw_read_count is the read count which you use as input for e.g. DESeq2. It has been estimated using RSEM
normalized_read_count is equivalent to the scaled_estimate from TCGA. This is the estimated fraction of transcripts made up by a given gene, as estimated by RSEM. Multiplying this value with 1e6 yields the TPM.


Is there any way of using biopython to write Swissprot files?

Using SeqIO.index rather than SeqIO.parse lets you read all the records into a dict, from which you can then extract the raw entry:

from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
for key in record_dict:
    print(record_dict.get_raw(key).decode())

Now you should be able to apply your test for a transmembrane protein to each entry, and write out only the ones you want to keep.

from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
out = open('Output.txt', 'wb')
for key in record_dict:
    record = record_dict[key]
    for i, feature in enumerate(record.features):
    if feature.type == "TRANSMEM":
        out.write(record_dict.get_raw(key).decode())


How to downsample some of the samples in RNA-seq data?

You have a few options:

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

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