Merge 2 VCFs from different variant callers

The merge utility from SURVIVOR seems to be what you are looking for. I did not try it yet, but it seems to be especially designed to handle SV.

The simplest thing to do would probably to first merge your patients files together for each tool and have one multi-patient file per variant caller. Then I would merge the multi-sample files of the different tools.

Here is a detailed blog post describing the different challenges associated with merging structural variants from different callers and explaining the steps they took to address it. All the code to reproduce the pipeline is available on a github repo.


Strategies to disentangle the effect of two genes on phenotype

Maybe you could fit a mixed model to predict phenotype, using the two genes as fixed effects and cancer type as a random effect, and ask the model for the separate effects of each gene, and the interaction between them.


Does rRNA depletion protocol give higher number of mapped reads in Intronic regions?

It's not so much that you have "intronic contamination" or "genomic contamination", rather you're not selecting explicitly for full-length mature transcripts with rRNA depletion. That is the most common cause for higher intronic read rates. There's nothing you can do about this post-hoc, just continue along.
BTW, many lncRNA's are polyadenylated, so you'll keep them with poly-A selection.


I think that there is a good argument for using technologies adapted from RNA-seq for metagenomic abundance quantification. The tool that I am most familiar with (that I use for this purpose) is kallisto, which uses pseudoalignment instead and is therefore less compute intensive.
For some background on the conceptual linkages between metagenomic abundance and transcriptomic abundance you can see the kallisto metagenome abundance paper or the explanatory blog post.
I think that it is preferable to use a tool that is specifically adapted for metagenomics, as there are always some numerical issues that might have been accounted for by the authors. I think that for this purpose it is maybe best not to naively use tools like mpileup, though I have no specific objection to it. However, if you look around you can see various papers using e.g. Cufflinks and Cuffdiff for effectively the same purpose, and Cufflinks is definitely intended for RNA-seq.
Depending on your specific application, you might want to use a different tool like mOTUS2, if you are interested in estimating abundance of specific genomes or lineages in a metagenome mixture. This method is probably more accurate as it is based on averaging across the whole genome, but it is not very useful if what you are interested in is more e.g. functional profiling of the metagenome.


From this document, it looks like REMARK 250 is the way to go.


log rank p value interpretation from gepia 2

The plot on the left shows that the gene has very different expression between two groups. However, according to the graph on the right, these high expression is not directly related to survival because both groups have similar survival time.
The gene expression so high might be to compensate for other deficiencies or problems.


Analyzing transcripts for non-coding RNA

ncRNAs are by itself not special. Sure they might be more unstable and less abundant than other transcripts but for the analysis it is more or less the same. You should take your transcript level counts and aggregate them to the gene level. The tximport package comes to mind. Every transcript (given it is annotated) is linked to a gene. ENST00000611544 is connected to the U2 gene but these snRNAs are also short and I doubt you see them in normal RNA-seq.
In the end you can process data as any RNA-seq. Quantify your data against a transcriptome, e.g. with salmon, then aggregate transcript level abundance estimates (so basically counts) to the gene level and see whether your genes of interest have counts you can work with or whether you mostly have zeros. A transcriptome can be obtaine e.g. from GENCODE. The available annotation files will tell you which genes these transcripts originate from. Check for the GTF files.
Towards enrichment analysis, well it all comes down to the annotations you have. If the ncRNAs are annotated in any pathway then you can perform standard analysis with tools such as goseq or gprofiler2. If they are not annotated, well then there is no information available.


I don't believe this statement is true:

To solve this issue, both TopHat2 and HISAT2 chop up the unmapped reads into around 25bp fragments.

TopHat2 did not support spliced alignment directly.  This is why it followed this general approach of chopping the read up into small parts, followed by "contiguous" alignment of the sub-reads.  Hisat2 / Hisat-genotype support spliced alignment directly insofar as they don't pre-partition the reads into small parts, but allow dynamic extensions of alignments across splice junctions (in a way similar to how STAR works).  This is generally a fundamentally better approach and is one reason, among many, why TopHat/TopHat2 have been put into maintenance-only mode and have been officially superseded by Hisat2/genotype.


Using column 2 of one file to match with two columns of another file, and append

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

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

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

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

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

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

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


How can I improve ESS values for ucldMean in beast?

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

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


Single Cell RNA seq Analysis for low input Cells?

You can perfectly use any of the single cell processing packages out there. Seurat, Monocle3, Scanpy, Kallisto, etc. Following some of the tutorials and vignettes will most likely help you figuring out how to proceed.

These packages provide useful functions to explore single cell data, process the data and each of them has a series of advanced features to address specific questions. You don't need 10k cells to look at number of gene/cell, transcript/cells, % mitochondrial genes, etc. You can perfectly normalize and standardize 500 cells and run PCA, and analyze the complexity that describes your data.

To address if 500 cells is enough for you pilot we need to know what question you want to address.
The following tool for example gives some intuition on how many cells to sample if you wish to detect n-cells per cluster.

https://satijalab.org/howmanycells

Hopes this helps


Increase number of threads for GATK 4.0 HaplotypeCaller

I am assuming you are using the non-SPARK version of the method. To specify the number of threads you wish to use with HaplotypeCaller, include --native-pair-hmm-threads (documentation). This will only parallelize the pair hidden Markov models (pair HMM) process. This is an argument for the method HaplotypeCaller itself, so it should be located after the method is called as opposed to one of the --java-options.

I do recommend trying out the SPARK implementations of the GATK tools. The whole method will take advantage of the threads you specify. I've been running SPARK locally (i.e. not on Firecloud) on computational servers. I've been impressed with the enhanced performance.


How to efficiently extract coordinates from a `pybedtools.Interval`?

The reason why unpacking an Interval results in strings instead of ints might be related to how indexing by position works.

Even without considering the string/int issue, it seems slightly more efficient to access the interval's attributes rather than unpacking them directly:

In [311]: %timeit [(chrom, start, end, strand) for (chrom, start, end, _, _, strand) in my_bedtool]
306 ms ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [312]: %timeit [(interval.chrom, interval.start, interval.end, interval.strand) for interval in my_bedtool]
238 ms ± 24.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

To resolve the ugliness issue, attribute access can be made more compact through the use of operator.attrgetter:

In [313]: from operator import attrgetter

In [314]: get_interval_attrs = attrgetter("chrom", "start", "stop", "strand")

In [315]: %timeit [get_interval_attrs(interval) for interval in my_bedtool]
239 ms ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


What is a samtools mpileup reference skip?

That documentation has been expanded recently (see PR #1055) and now describes the characters seen in this column in more detail:

Forward   Reverse                    Meaning
---------------------------------------------------------------
 . dot    , comma   Base matches the reference base
 ACGTN     acgtn    Base is a mismatch to the reference base
   >         <      Reference skip (due to CIGAR "N")
   *         *      Deletion of the reference base (CIGAR "D")

(And so on — see the pull request or doc/samtools-mpileup.1 on the develop branch for the complete pileup format description. The documentation on that web page will be updated when this is released in a new samtools release.)

As noted in https://bioinformatics.stackexchange.com/a/421/134, SAM has always had separate D (deletion from the reference) and N (skipped region from the reference) operators. These are treated similarly computationally but emitted intentionally by aligners to distinguish introns and the like from other deletions.


ERROR: (subscript) logical subscript too long

I am answering my own question here.

I realised I should have been using

rawCounts <- rawCounts[rownames(rawCounts) != geneID, ]

Problem solved.


After artificially creating events in a FASTA file, how do I keep track of the old coordinates?

I've written a handful of programs from scratch to simulate mutations and variations in real or simulated sequences.

The trick has always been to sort the variants by genomic coordinate, apply the variant with the largest coordinate first, then apply the variant with the second largest coordinate, all the way down to the variant with the smallest coordinate.

Indels and other variants that affect sequence structure only affect subsequent coordinates. So going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence.


ComBat for batch effects removal on scRNA-seq data

Look at this recent paper that uses ComBat on scRNA-seq data for batch effect removal and states that it "successfully does so".

I also suggest that you check out this publication on Distribution Matching Residual-Nets. Authors evaluated their method also on scRNA-seq data and thus it may be something you are looking for. I personally played a bit with their code and found it straightforward to use.


How to Convert From Gene Location to SNP Location

You first need to get the coordinates of the "cytoband" 21q22.11, which I would search for using the UCSC table browser. You need the following settings:
Clade: Mammal
genome: Human
assembly:
group: Mapping and Sequencing
track: Chromosome Band
table: cytoBand
region: position chr21
output format: BED

(you can leave other options as default)

That will give you a bed file of cytobands on chr21. You'll need to search (e.g. grep) for your cytoband: q22.11.


How to manipulate a reference FASTA or bam to include variants from a VCF?

GATK has a solution that might work for you:
FastaAlternateReferenceMaker, which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records."

  Input

  The reference, requested intervals, and any number of variant ROD files.

  Output

  A FASTA file representing the requested intervals.

  Usage example

java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -L input.intervals \
   -V input.vcf \


Each color is a clone (set of mutations) or an individual mutation (depends on how the input data was generated). This is a cancer with the bottleneck likely being induced by a treatment (e.g. UV, chemo, drugs, etc...)

Based on the vignette this is likely an AML (leukemia).

It looks like grey are healthy cells. At day 0 a specific clone (red) existed at a low frequency (~20-40%). The patient (or cells) was then exposed to a treatment and on day 150 another measurement was taken. In this case two novel clones (subclones of red) have taken over the population. They seem to co-exist or are co-dominate clones with orange being 50-60% of the population and yellow being 30-40%.

The example says "Simple plot with multiple subclones" so: the orange and yellow populations have additional (likely mutually exclusive) mutations that were not present in the red clone - but biologically also have all the mutations in the red clone. So my first assumption was right - each color is a proxy for a set of mutations.

You could have generated the mutation data (i.e. allele frequencies) from a wide variety of sources (scRNA-Seq, exome-sequencing, PCR genotyping, etc...)


How to retrieve the DNA sequence from a particular species and region through command line?

Samtools
Samtools is a toolkit that includes a binary called samtoools, which can be used to query a FASTA file — to extract a subsequence based on some coordinate range you specify.
If you do not have samtools installed, you need to add it to your computer. For OS X, you could use Homebrew, for example:
$ brew install samtools

For Ubuntu, you could use apt-get:
$ sudo apt-get install samtools

Reference assembly
The next step is to download a reference assembly for your genome and genomic build-of-interest.
There is a research group dedicated to making Xenopus datasets available. If you were interested in v10 of the Xentr genome, for instance, you could do the following or similar to download it:
$ curl --output "XENTR_10.0_genome.fasta.gz" "http://ftp.xenbase.org/pub/Genomics/JGI/Xentr10.0/XENTR_10.0_genome.fasta.gz"

You could also dig through the NCBI site to get a similar link. It's the same idea — just change the web address, accordingly.
Indexing
Indexing makes it possible to extract an arbitrary subsequence of interest from this FASTA file. You need to index, but you only need to do it once:
$ samtools faidx XENTR_10.0_genome.fasta.gz

This creates two new files: XENTR_10.0_genome.fasta.gz.fai and XENTR_10.0_genome.fasta.gz.gzi.
Note: You need to keep these two files in the same directory as the original, compressed FASTA file XENTR_10.0_genome.fasta.gz.
Querying
Now you are ready to query, or extract a subsequence of interest. If you wanted the first hundred bases of the first chromosome, for instance:
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100
>Chr1:1-100
cctcagcgtcttcttctttctgctcccactgcattagcatagggagagagggccacagag
gctaggtagcattgccggctagcaatcttcagcgtcctct

The search term used here is Chr1:1-100. Replace Chr1 with the chromosome of interest and 1 and 100 with the start and stop coordinates, respectively, of interest. Coordinates used with samtools queries use a one-based indexing scheme.
Note: This uses the chromosome naming scheme in the XENTR_10.0_genome.fasta.gz file (Chr1, etc.). Other research groups (NCBI) might use different chromosome names. If you get your FASTA from somewhere else, investigate the FASTA file to determine the naming scheme being used.
If you want to redirect output to a file, use the output redirect operator (>):
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100 > mySubsequence.fa


A member of the SRA submission staff pointed out that using

prefetch --type all SRR5253957

will download the original files. In this case, it means running the above within an EC2 instance colocated with the S3 bucket (so, us-east-1) and having installed and configured SRA Toolkit to work from AWS (as per this documentation).

Unfortunately, the particular files I am concerned with are not currently accessible, but this should generically work in most situations.


Error in annotating row names as data points in a PCA plot

Maybe try

plot(evecdat$PC1, evecdat$PC2, xlab="PC1", ylab="PC2", pch="24", cex.lab=2.5, cex.axis=1, type="n")
text(evecdat$PC1, evecdat$PC2, evecdat$Sample, col=as.vector(evecdat$color))

This should plot the sample names as text (taken from the evecdat$Sample column), colored as specified in the evecdat$color column. See ?text for details.


How to set the position of groups in a Seurat object on a FeatureHeatmap plot

I don't think this is possible in Seurat v2, but in v3 you can change the factor levels of the grouping variable to change the plot order:

library(Seurat)

FeaturePlot(object = pbmc_small,
            features = head(VariableFeatures(pbmc_small), 2),
            split.by = 'groups')

Change the order:

pbmc_small$groups <- factor(pbmc_small$groups, levels = c('g2', 'g1'))

FeaturePlot(object = pbmc_small,
            features = head(VariableFeatures(pbmc_small), 2),
            split.by = 'groups')


Normalize ATAC-seq/Dnase-seq sequencing reads coverage signals over estimated background

Supplementary File 5 and 6 contain the code you're looking for.


How to specify input resfile for Rosetta?

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

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

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

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


how do I predict bacterial small non-coding RNA for a specific mRNA?

Answer from @llrs, converted from comment:
You could try to download RNAseq data from this bacteria and see if some of the RNA align to your sequence. You could try to blast your gene and look for the sequences that have low identity similarity (or classify them by type and then look what is each hit).


How to quickly determine mutations in a read of a sam file?

For qualitative analysis, you're probably better off using something less granular like IGV or IGB. However, if you really want to look at a couple of reads:

If you're willing to ignore sequencing errors, you can inspect the CIGAR string or the MD tag, both of which give information on the alignment of a single read.

The CIGAR string gives details on insertions, deletions, clippings, matches and mismatches. From Genome Analysis Wiki,

  The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference. For example:

RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

Most common usage of the CIGAR string uses M (match/mismatch), I (insertion), D (deletion), S (soft clipping), and H (hard clipping). Note that = (match) and X (mismatch) are available as alternatives to the less informative M, but they are less widely used.

The MD tag gives specific details on mismatches and deletions. From the SAMtools tags specification,

  The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string ‘10A5^AC6’ means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which is different from the aligned read base; the next 5 reference bases are matches followed by a 2bp deletion from the reference; the deleted sequence is AC; the last 6 bases are matches. The MD field ought to match the CIGAR string.

Note that neither of these will give you any idea of structural variants in short reads, and neither will be particularly readable (or helpful, due to the higher error rate) in long reads.


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

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

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

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

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


Fast5 files need to be converted to fastq first prior to mapping. There are currently no tools available to map the entire raw signal directly to a reference, although a few tools are available for identifying short sequences within the raw signal.


How to make a video of a gromacs .trj file with VMD

No need to use gmx dump, you can simply load a xtc or trr trajectory in VMD using:
vmd path/to/gro_file.gro other/path/to/trajectory_file.xtc

As always you can also use the GUI to first open your gro file, then right click on the molecule and select Load data into molecule before selecting your xtc or trr file.
After putting the representations as you see fit, if you want to export this as a movie (GIF, MP4, others...) you can go in Extensions -> Visualization -> Movie Maker. See Trajectories and Movie Making VMD official tutorial for more details.


Which of these enzymes is a better target for virtual screening?

The ultimate choice stems from your individual interests and not a scientific reason... but having said that some pointer of issues may be handy.

I am assuming that you want to make an inhibitor that messes with the activity and not a ligand that binds somewhere and you then plan to add a moiety to make it a protac to get it degraded or to make it fluorescent to mark it. If not, ignore any talk of activity.

Crystal

is a crystal structure available?
does the structure have the cofactor if any?
is a structure available (itself or homologue) with the native metabolite bound?
how is the active site? In the case of many cytochrome P450 enzyme the entranceway is really deep and is guarded by a flappy gate, both of which confer a great deal of specificity regardless of the residues in the active site, which are often just water coordinating. Your first enzyme, is that a NADPH/FAD or a cytochrome enzyme?
What is the mechanism and is there a classic type of inhibitor chemistry? I.e. catalytic triad with a cysteine nucleophile is a nice target.(1)

Methodology

Wet lab: A problem with computational chemistry/biochemistry is that someone needs to make the protein and physically test your compounds. Are these easy to express? This may sound like not your problem, but a lot of issues can come from this latter stage.
is there a close homologue you want to avoid targetting? If so, you'll have to do twice the docking (not a big deal and potentially more interesting). (2)
is the method you plan to use okay with covalently bound ligands? Dopa decarboxylase sounds like a PLP enzyme. If it is, making an inhibitor using classic functional groups that nucleophilically attack the pyridine ring is very straighforward, but not all docking methods are okay with covalent ligands.
Are any hits likely to follow Lipinski's rule of five? If you have a very charged active site, your hits will be metabolite-like (polar) and not drug-like (horrendously insoluble).

A minor thing is to make sure that nobody else is working on it or you may be left behind. Industry is all cloak and dagger, but the academic SGC has a list of [target enabling packages)[https://www.thesgc.org/tep] and their future ones (nominations).

Physiology

This was most likely in your lit search. But I am double checking.

You have chosen these enzymes because their unchecked activity is disease-causing, right?
And they do not have any other roles in other tissues?
Abolition of activity does what you want?

Other hits

Lastly, but not least. See which of these have many paper about barely active hits: this should mean it is hard. For example, bacterial alanine racemase has been a holy grail for a while, but progress is slow as it has been described as a "ship in a bottle" due to its long and unforgiving entranceway.

Addenda

A catalytic triad is acid-base-nucleophile: if you have an electrophile with a leaving group (including internal ones, such as epoxydes) you can virtually screen more complex small molecules that resemble your ligand, but with the difference that there is a leaving group, resulting in an irreversible inhibition. Whereas competitive inhibitors require nanomolar affinity to outcompete the natural ligand, an irreversible inhibition can even slightly worse and will work.
The problem of homologues.

Paralogues (homologous genes in the same organism): you want to target protein 1, but you do not want to target protein 2. This is a common problem with kinases. You therefore need to find a hit that works better in one than in the other, which with a big of rational inference is easy.
Homologues (homologous genes in different organism). You want to kill species A (say a fly), but not B (say a human). Same thing.

Twice the docking. It is more of a challenge, but makes for a better spin on the story. There is an explosion of people doing virtual screens that do not go anyway afterwards, so it is getting harder and harder to publish studies without extensive empirical backing. Hence the need for an interesting angle.


Percentage distribution of cells in all clusters based on their treatment condition?

Here is a solution using dplyr and ggplot2:
library(Seurat)
library(dplyr)
library(ggplot2)

meta.data <- pbmc_small[[]]

# create random classifications for the sake of this example
meta.data$condition <- sample(c('A', 'B', 'C'), nrow(meta.data), replace = TRUE)

counts <- group_by(meta.data, condition, res.1) %>% summarise(count = n())

ggplot(counts, aes(res.1, count, fill = condition)) +
  geom_bar(stat = 'identity')


As explained in the pep docs, the sample_name column should be considered a unique identifier for a sample. The PEP requires that One row corresponds to one sample You have multiple rows with the same name; this is not allowed.

You simply need to recast your annotation sheet such that it has file1 and file2 columns, but only 2 rows.


What does liability mean in GWAS heritability?

When you use linear mixed models to estimate heritability you assume that the underlying trait is normally distributed which is called a disease liability scale.

For continuous traits this is not a problem but for binary traits, this becomes an issue because you have a 0/1 value for a phenotype and usually there is a higher proportion of cases in the study sample than the general population prevalence of disease, which leads to an ascertainment bias.

The population prevalence also varies between populations, for example the prevalence of malaria in one continent is vastly different from another.

So to make the heritability estimates comparable and also to correct for the ascertainment bias, the observed scale heritabilities for dichotomous traits are transformed to the liability scale taking population prevalence and the case proportion into account.

Some good papers for reading:

Falconer D.S. The inheritance of liability to certain diseases, estimated from the incidence among relatives. Ann. Hum. Genet. 1965;29:51–71
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3059431/


What is the meaning of these misaligned reads in a sequencing run?

The unaligned parts of the bottom two sequences are from elsewhere in the sars-cov-2 genome if you blast them. Likely causes of this are library prep. artifacts (not terribly common, but you'll find a few cases here in there once you have millions of reads) where fragments end up ligating together or (probably more likely) rearrangements of the viral sequence. My guess would be that the latter is what occurred, given how rapidly viruses reproduce. Whether that arrangement can lead to viable viruses I don't know (that's too far out of my field).

Assemblies represent at best a snapshot of a single example of a species at a single moment in time. Actual individuals in the wild can vary quite substantially, especially when they don't need to be viable.


How to fragment genomes into non-overlapping sequences of differing sizes?

You could use numpy.random.randint to sample integers from a uniform distribution, from 1000 to 15000 nt, extracting subsequences of that length from a linearized multi-FASTA file.
Here's some code:
#!/usr/bin/env python

import sys
import numpy as np

low = int(sys.argv[1])
high = int(sys.argv[2])

def process_record(header, sequence, low, high):
    next = 0
    fragment = 1
    while next < len(sequence):
        new_next = np.random.randint(low, high) + next
        if new_next > len(sequence) - high:
            new_next = len(sequence)
        sys.stdout.write('{}:{}\n{}\n'.format(header, fragment, sequence[next:new_next]))
        next = new_next
        fragment += 1

header = ""
sequence = ""
for line in sys.stdin:
    if line.startswith('>'):
        if len(sequence) > 0:
            process_record(header, sequence, low, high)
            sequence = ""
        header = line.rstrip()
    else:
        sequence += line.rstrip()

if len(sequence) > 0:
    process_record(header, sequence, low, high)

To use it:
$ ./13879.py 1000 15000 < in.fa > out.fa

Fragments will be disjoint (non-overlapping, contiguous) subsequences.
I add a check to determine if the last fragment of a FASTA record would be shorter than low (1kb) in length, and, if so, append that fragment to the final subsequence. It is possible for that last fragment to be up to 16 kb - 1 nt in length, given your parameters:
if new_next > len(sequence) - high:
    new_next = len(sequence)

You could change this logic to discard the last fragment (i.e. break), or simply allow fragments less than 1kb in size, depending on what you need.
EDIT
Here's a slight modification that should allow to you preserve fragment sizes within your bounds.
A cartoon may help illustrate the problem. Let's set up a FASTA record, which has been broken up into n fragments:
   1       2           n-2      n-1      n
|-----|---------|...|-------|---------|-----|

We may get unlucky and sample a random integer for fragment n-1 (frag(n-1)) which makes the length of frag(n) less than low (1 kb in your example).
In my first suggested solution, frag(n-1) and frag(n) are appended to one another. This is why you get some fragments longer than high (15 kb).
Another option is to discard frag(n), which is of length less than low (1 kb). Here's some code:
#!/usr/bin/env python

import sys
import numpy as np
import timeit

low = int(sys.argv[1])
high = int(sys.argv[2])
in_fn = sys.argv[3]

"""
Globals
"""
fragment_lengths = []

"""
Process a record
"""
def process_record(header, sequence, low, high):
    next = 0
    last = 0
    fragment = 1
    fragment_sequence = ""
    while last < len(sequence) - high:
        next = np.random.randint(low, high + 1) + last
        fragment_sequence += sequence[last:next]
        # Comment the sys.stdout.write and uncomment the fragment_lengths line, if measuring fragment length distributions
        sys.stdout.write('{}:{}-{}:{}\n{}\n'.format(header, last, next, fragment, sequence[last:next]))
        #fragment_lengths.append(next - last)
        last = next
        fragment += 1
    """
    Uncomment to compare input sequence with fragments
    """
    #sys.stderr.write('In > {}\nOut > {}\n'.format(sequence, fragment_sequence))

"""
Parse records
"""
def parse_records(in_fn, low, high):
    header = ""
    sequence = ""
    with open(in_fn, 'r') as ifh:
        for line in ifh:
            if line.startswith('>'):
                if len(sequence) > 0:
                    process_record(header, sequence, low, high)
                    sequence = ""
                header = line.rstrip()
            else:
                sequence += line.rstrip()
        if len(sequence) > 0:
            process_record(header, sequence, low, high)

def count_elements(seq):
    hist = {}
    for i in seq:
        hist[i] = hist.get(i, 0) + 1
    return hist

def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped

def main(in_fn, low, high):
    """
    Comment below if testing the distribution of fragments
    """
    parse_records(in_fn, low, high)
    """
    Uncomment to test the distribution of fragments from
    running fns on same input 100k times
    """
    #wrapped = wrapper(parse_records, in_fn, low, high)
    #timeit.timeit(wrapped, number=100000)
    #print(count_elements(fragment_lengths))

if __name__ == "__main__":
    main(in_fn, low, high)

If frag(n) is needed, and you need to preserve the uniform distribution (or whatever fragment length distribution you want to model), then rejection sampling might be an option.
With rejection sampling, you keep fragmenting your input sequence, until one of the distributions of lengths gets you a set of fragments that entirely cover the input sequence.
If you need that, follow up and I'll take another look.


I've had a quick look at the file you linked to and this is a featureCounts file - this list the number of reads for each gene/transcript (COUTNS.txt in the file name pretty much gave it away that it was going to be a file of, well counts).
I can't help you with Galaxy, but if you want to do an RNA-seq analysis, the count data + metadata is likely all you'll need.
If you really want to do the preprocessing and mapping yourself, the fastq files can be downloaded from the SRA archive (SRP242087) - click on the SRA link in your whole series link.


This could be organism-specific. We don't have a lot of info so far, so I would check a few more things:
Run something like FRC_align. Check if there's a clear signal between regions flagged as suspicious by it and your coverage graph.
Is it a eukaryote? Plant? Check where mitchondria and chloroplasts are on the plot. They will have different GC/coverage signals than the rest of your assembly, which is fine, expected even.
Worried about contamination? Run kraken2 on your raw data and see what you get there.
Run RepeatMasker and check if repetive regions contribute to the effect you see.
Finally, unless your specimen is low heterozygosity, haploid, or inbred, you probably have some effect from ploidy. There are numerous ways to check that. A self-self dot plot with MUMmer and looking coverage histograms are a good start. The idea is that you should see two peaks (let's assume diploid). One for the split haplotigs, and one for collapsed.
I checked the GC distribution of 45 genomes in your order. It's kind of all over the place, but there are some bimodal ones:

What does GC vs Length look like? Here's the one from the 45 genomes mentioned above:


blasting a refseq protein does not show the protein in the result set

NP_420767 is represented by the non-redundant refSeq protein WP_010919826, which has the same amino acid sequence. This is not very clearly annotated, but if you scroll down to the sequence in the GenPept entry for NP_420767, you'll see the following:

CONTIG      join(WP_010919826.1:1..799)

See here for more info.


How to confirm exon shuffling in a gene?

These suggestions are based on your comment that you have the gene sequence.

For a statistical approach, you could formulate a null model and show that the sequence similarity you observe is much higher than expected under the null model, but reasonable under your exon shuffling model. Ideally, you would show this individually for each exon, to convince the reader about the exons' disparate origins. I'll propose a reasonable null model, though suggestions are gladly accepted here. Pick a random gene, assume it's built using exon shuffling, and search the rest of the transcriptome for sequences that best substantiate that argument.

For an evolutionary approach, you could search for orthologs (homologues? I never remember the difference) of your gene in other species. Hopefully, the results would show that your gene formed more recently than the exons it's built out of. If your gene actually has some other origin, maybe you'll find out it's older than some of its putative components. This can't exactly "confirm" your hypothesis, but it will test your hypothesis and just be generally a cool thing to try.


Where can I get ensembl Medaka genome for RNAseq

Download the ones without a chromosome number in the filename. Those files are the whole genome.


I would strongly recommend to use Bioconda instead of compiling. It is one of the channels ("distributors of packages") in Anaconda Cloud which is basically a package source for conda.

In short - conda allows you to create virtual environments isolated from system and to install various precompiled software packages with one command (so, by the way, you can easily use multiple software versions on one system). Also, Bioconda is just one channel aiming on various bioinformatics software, but many more packages can be found in aforementioned Anaconda Cloud.

Tutorial how to start using Bioconda is on its homepage, and it is really easy (far more easier than to struggle with compiling).

Bioconda HMMER package

P.S. Of course there are cases in which you need a specific build of given software. Than it is really needed to compile manually.


Converting VCF file to PLINK bed/bim/fam files

Given that Plink reads in VCF files natively, and Plink is your desired output format, I expect that plink --vcf <file> would be the best option.

There are caveats associated with the types of variants that Plink will pull in from the VCF files, mostly related to limitations of Plink's own file formats. But if your end-goal is Plink, then that should be fine.


This was answered by the OP -- it was a bug that has been fixed in a subsequent release, 4.1.5.0. See the release notes here: https://github.com/broadinstitute/gatk/releases/tag/4.1.5.0

  Made callable sites a Long to avoid integer overflow (#6303)

Issue ticket: https://github.com/broadinstitute/gatk/issues/6302


Number of reference sequences in a SAM file

There are multiple ways. Here are two:

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

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

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

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

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

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

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

Or:

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


How to run Cox proportional Hazard model for each group in R?

The fact that you supplied GROUP as a covariate, you are getting 9 coefficients as one of the ten levels is used as the reference level.

If you would like to fit one model per group (and in this scenario you would not be correcting for GROUP), you should first split your data frame and fit your models on the resulting list's elements:

lapply(split(DF, DF$GROUP),
       function(x) coxph(serve_object ~ age + Smoking + Diabetes, data = DF))

Example with the mtcars data set when grouped by number of cylinders, cyl having values of 4, 6 or 8:

> lapply(split(mtcars, mtcars$cyl),
+        function(x) lm(disp ~ hp, data = x))
$`4`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
    59.0369       0.5579

$`6`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   290.9178      -0.8799

$`8`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   320.2083       0.1572


Description of how to use biomaRt here. Filter by your list of gene IDs. Get the descriptions etc as attributes. Eg:

> library(biomaRt)
> ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
> IDs <- c("BRCA2","BRAF")
genedesc <- getBM(attributes=c('external_gene_name','description'), filters = 'external_gene_name', values = IDs, mart =ensembl)
> genedesc
  external_gene_name                                                                      description
1               BRAF B-Raf proto-oncogene, serine/threonine kinase [Source:HGNC Symbol;Acc:HGNC:1097]
2              BRCA2                  BRCA2, DNA repair associated [Source:HGNC Symbol;Acc:HGNC:1101]


What does the 'suppressed' folder mean on the NCBI ftp genomes website?

Here you have the explanation:

Suppressed genomes: Several hundred RefSeq bacterial genomes that did
not pass assembly or annotation quality validation have been
suppressed(RefSeq genome quality criteria). An FTP report file listing
the most recent suppressions will be provided soon in the RefSeq FTP
site.


Import gene list to Seurat to define cell types

Seurat does not define cell types by name. It clusters and assigns each cell to a cluster, from 0 to X.
If your data has the cell type (e.g. B,T, Mast cells) it means that someone annotate the clusters so that they have a biological meaning.

You can assign different names to the clusters by using the AddMetaData function.

sampleID  cellType
cell1     B
cell2     B
cell3     T
cell4     T

newMeta <- oldMEta
newMeta$origType <- newMeta$cellType
newMeta$cellType <- c("B1","B","Treg","Tkill")

object <- AddMetaData(object, newMeta)

The new metadata should look like

sampleID  origType  cellType
cell1     B         B1
cell2     B         B
cell3     T         Treg
cell4     T         Tkill

You can then swap the labels using SetAllIdent(object = test, id = "cellType"). Your tSNE or uMAp, will now show the new cell types.

Have a look at Manually define clusters in Seurat and determine marker genes


You have to perform an alignment first. Phangorn use phylogenetic estimation and this requires a precise alignment position to track homologous nucleotide mutation. Maximum likelihood theory will assess the probability of a given tree for each homologous nucleotide site, if they are not alignment it can't do it and will violate the assumptions of phylogeny.

If the sequences are different lengths, which it means both nucleotide and indel ('-') sites, it means there hasn't been or there is an error with the alignment. I use muscle these days having previously used clustalo. Its really easy to use and easily configures on a linux server and its super quick. If you get stuck let me know.

I had a look at your alignment and can verify it hasn't been aligned.


Generating DNA sequences with constraints

You can use a $l$-order Markov chain. Here is the procedure:

Count $l$-mers in your genome. For small genomes, you can do that in Python. For large genomes, you may need jellyfish or KMC3.
Draw a $l$-mer randomly based on the distribution of $l$-mers. This generates the first $l$ bases.
Let $s$ be the last $(l-1)$ subsequence from the generated sequence. There will be at most four $l$-mers prefixed with $s$. It gives you a distribution of the next A/C/G/T base conditional on $s$. Randomly draw a base from this distribution. You get the next base.
Repeat step 3 until you generate a sequence long enough.

The generated sequence will have the same $l$-mer distribution to the original sequence. The GC content will be the same, too. In your case, if $k$ is small enough, you can let $l=k$. If $k$ is large, you have to choose a smaller $l$.

The procedure above doesn't consider the strand symmetry of DNA, but probably this doesn't matter much. It is also possible to adapt the algorithm to address this.

As to implementation, you can use python, but a high performance language like C will be tens of times faster at least.


How to install genometools with anaconda?

genometools-genometools seem to be the current version. By installing the genometools package I received the same error like you. I guess it's because the binary isn't copy to bin during installation.

genometools-genometools seems to work fine. And it has 2k (2000) and not 2 :)


Disclaimer: I'm a developer for http://www.sequin.xyz.

Sequin is a new set of spiked-in controls for next-generation sequencing, and that includes Illumina. We design controls for RNA-Seq, genome sequencing, metagenomics etc.

Please study our papers if you want more details.

Reference Standards For Next-Generation Sequencing by Hardwick should be a good start.

Here's my list:

Measured counts against input concentration. Pearson's correlation, spearman correlation and regression slope.
PCA analysis
Whether the synthetic aligned to the synthetic reference
Spiked-in amount (e.g. dilution)
Diagnostic statistics
ROC Curve
Detection limit

(1) is important because the controls should give you correlation = 1, and slope = 1 in a perfect flawless experiment. This is a sample output (you can read more about it in our papers):

(2) This is our sample PCA (our paper has the details)

(3) is also important but you can only do that in simulations.

(5)

(6)

(7)


How do I interpret my phylogenetic tree construction from Glucosidases in the Actinobacteria class?

Bootstrap values usually give you useful information, be they low or high. Usually, a "high" bootstrap value is something >= 70%. That means a clade is strongly supported; lower values are telling you the clade is not clearly separated from others. This could be for several reasons, e.g.:

1) A subpopulation is not really separated from others, e.g. the viral sequences in  town "A" and town "B" are completely admixed, meaning the population is panmictic and there is continuous gene flow between the two of them;

2) the populations diverged not long ago, it's too soon to have any phylogenetic signal, at least using the marker you are using (different genes evolve at different rates). Usually, housekeeping genes (which are usually very conserved) are used for populations that supposedly diverged long time ago.

The interpretation of your tree should take into account the biology and epidemiology of your target organism. If you really want to learn more, I suggest a lot of reading, and definitely a class dedicated to molecular phylogenetics.

This might be a good start:
https://www.amazon.com/Phylogenetic-Handbook-Practical-Approach-Hypothesis/dp/0521730716

Also, the NJ algorithm should allow for uneven branch lengths, maybe the tree you are showing is the consensus tree?

Hope it helps

Max


Picking number of samples for clustering

In case there is no systematic difference between the samples with or without the clinical data (i.e. one set is microarray and the other is RNA-seq), I would use all of the samples. And if the cluster labels obtained with "clustering 50 samples" and those obtained with "clustering 100 samples" does not differ too much, then doing either of two would not have much effect anyway. If you would like to compare cluster labels in that way and assuming you would want to test a number of ks, it would be better to compare each combination of ks and number of samples.


Understanding ViennaRNA RNAdistance scoring table

Ok, the answer was right here and I missed it, posting it in case anyone else will do the same mistake as me in the future.
U = unpaired
P = paired
H = hairpin
B = bulge
I = internal loop
M = multiloop
S = stem (or stack)
E = external loop


Sibling vs Parent/Child relationship detection using whole-genome sequence data

With at least three samples, you can determine this by looking at the IBS0 output column (from either plink2 --make-king-table or KING): if the relationship in question has a very low IBS0 value compared to the other pairs, that's parent-child.  You can try to do the same thing with only two samples, but with no other IBS0 values to compare against it's easier to misjudge.


well identifiers in a 1536-well mtp plate

and here's a scheme I made up myself:

AA01..AA48
..
AZ01..AZ48
BA01..BA48
..
BF01..BF48

Pros:

alphabetic sorting of well ids does what you expect
it's clear what's the row and what's the column
it's easy to distinguish well ids of 384-well plate from those of a 1536-well plate

Cons:

it's not used anywhere AFAIK


correlation coefficient versus DEGs analysis: what's the best approach for low expressed genes?

First of all, RNASeq is extremely sensitive to batch effects.  Are these matched controls processed at the same time as the tumors?  IF they are not, lots of your differences will be caused by batch effect, and not cancer.

I believe the best strategy is to merge the FPKM columns and perform
correlation analyses of my gene of interest (GOI) versus all other
genes in both cancer and normal tissues.

What you propose isn't publishable.  Use the software that people have been using and validating for years.  Use either DESeq2 or EdgeR, or limma if you don't have access to raw gene counts.
If the underlying counts are low, the correlation won't mean much, because your real values could be quite off. If your gene has 2 counts in one sample, and 3 counts in another, that does not mean it really correlates with a gene that has 500 reads in the one sample and 750 in the other. I don't know that you can trick your way around it; you need more reads, or more samples for more power.


How to create a .phy file for constructing a phylogenetic tree?

You simply use Biopython's SeqIO and AlignIO.
The documentation is https://biopython.org/wiki/AlignIO
etc ...
The format you seek is phylip


Range overlap python error with genomic regions

You're reinventing bedtools intersect (or bedops), for which there's already a convenient python module:

from pybedtools import BedTool

s3 = BedTool('s3.bed')
s4 = BedTool('s4.bed')

print(s4.intersect(s3, wa=True, wb=True, F=1))

The wb=True is equivalent to -wb with bedtools intersect on the command line. Similarly, F=1 is the same as -F 1.


How to analyze co-occurrence of multiple SNPs?

I found the LDheatmap package in R very helpful at calculating linkage disequilibrium. You can make plots such as this:


Looking for a tool to find 16S rRNA in hundreds of genomes

You can use Barrnap

barrnap -k bac --threads "$(nproc)" \
        -o output_rrna.fna < input_file.fna \
        > output_rrna.gff 2> /dev/null

Put the file in one line of sequence per fasta id and get sequence with 16S annotation.

perl -pe 'if(/\>/){s/\n/\t/}; s/\n//; s/\>/\n\>/' output_rrna.fa | \
    perl -pe 's/\t/\n/' | tail -n+2 | cut -d\: -f1 | \
    grep -A1 16S > output_16S.fasta

Check the output_rrna.gff and make a BLAST of the sequence to be sure about the predicted 16S.


Is there a resource where I can donwload E. coli assemblies grouped by their pathotypes?

If ENA archive contains such information, then you can access it programatically through the EBI api.

Also you can access the NCBI database programatically through their eutils.

In both cases you can access annotations, though not sure if they contain what you are looking for in all the cases.


Is there a way to retrieve ENSEMBL IDs from a search query?

I would download the complete geneset FASTA file from ensembl and then parse out the genes you want with a script.

Alternatively, you could do similarly with API. Get all the gene names, select the ones you want, and then pull down each sequence.


Generating list of known fusions for Arriba from COSMIC

I managed to solve this with a one-liner:

grep -v 'Sample ID' CosmicFusionExport.tsv | cut -f 12 | grep -v -e '^$' | perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' | sort -u > CosmicFusionList.tsv

The various parts of the piped command are as follows:

Remove header:grep -v 'Sample ID' CosmicFusionExport.tsv |
Extract 12th field (Translocation Name):cut -f 12 |
Remove blank entries:grep -v -e '^$' |
Extract gene symbols from 12th column:perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' |
Remove duplicate lines:sort -u > CosmicFusionList.tsv
Optionally the final two stages could be:sort -u | gzip -c > CosmicFusionList.tsv.gz to make use of file compression (for what it's worth on a 3.9 kB file)

Since the authors of arriba have neglected to document this necessary data munging step, I've posted it here in Q and A format in the hope that it will be useful to others, and other solutions can be proposed.


Calling GEOquerry error while loading the library

rlang 0.4.10 is the current version at CRAN: https://cran.r-project.org/web/packages/rlang/index.html
Try to install it via install.packages, not via GitHub directly.


I have noticed a few errors:

1) You define a bunch of variables (Species, Class, ...) but then instead of creating your data frame data with these you create it by reading from text. I don't think this is good practice, you should do so like in your example code.

2) With the way you define, your data frame does not contain columns like Class or Species but columns V1 to V5, you should use argument header = TRUE in your read.table() call.

3) Even when you read your data as in step 2, you would need to transform your data table into "long format", please see below.

4) In your ggplot() call, you assign 7 to your y aesthetic, this does not make sense. You should either provide a vector of length equal to your other aesthetics (x and fill) or better a column name of values of interest.

# This is your code minus the redundant variables that you have created but are
# not used

library(tidyr)
library(ggplot2)

# Pay attention to the header argument
data <- read.table(
  text = "Class Cames1 Came2 Cams3 cams4
     EVV 2194 5435 2706 2902
     MaLR 503 108 184 231
     psy 284 187 413 295
     Nargo 274 1977 413 2195
     Pao 66 29 78 97
     Unknown 286 283 273 377
     Undefined 76 77 80 95",
  header = TRUE
)

# data in wide format
> data
      Class Cames1 Came2 Cams3 cams4
1       EVV   2194  5435  2706  2902
2      MaLR    503   108   184   231
3       psy    284   187   413   295
4     Nargo    274  1977   413  2195
5       Pao     66    29    78    97
6   Unknown    286   283   273   377
7 Undefined     76    77    80    95

# Making the data wider: Keeping Class column as it is but condensing the data
# within the rest of the columns into two newly created columns, Species_ and
# Value_

data_long <- gather(data,
                    key = "Species_",
                    value = "Value_",
                    -Class)

# First few rows of your data in long format
> head(data_long)
    Class Species_ Value_
1     EVV   Cames1   2194
2    MaLR   Cames1    503
3     psy   Cames1    284
4   Nargo   Cames1    274
5     Pao   Cames1     66
6 Unknown   Cames1    286

# x and y aesthetics are the newly created columns Species_ and Value_
ggplot(data_long, aes(fill=Class, y=Value_, x=Species_)) + geom_bar(position="dodge", stat="identity")


What are "split reads" and "intron clusters?"

split reads - These are read that have two or more alignments to the reference from unique region of the read.  In this example a 150bp read sequenced from RNA could have base 1-75 aligning to the 3' end of exon2 and bases 76-150 aligning to the 5' end of exon3.  This would be a split read because it have two alignments (exon2 and exon3) and those alignment are from unique (non overlapping) parts of the read, 1-75 and 76-150.

Chimeric reads are an example of split reads.  Many people use "chimeric reads" to describe a type of noise that occur during PCR or sequencing where the bases from two different molecules are reported as one read and the alignments of that read become split when it aligns back to the location of the two original molecules.

"intron clusters" - In this case I think that the LeafCutter software is describing stacks of split reads to cluster to support the detection of an intron.


Is it possible to filter contaminated reads for raw PacBio sequences (not HiFi reads) before assembly?

I don't think that in principle the issue is any different for PacBio and Illumina; it's just that the process might be noisier with higher error rates. Note that PacBio errors aren't randomly distributed substitutions and probably many (most?) short-ish k-mers (<30nt) will likely be intact. I would suggest doing k-mer filtering of your reads as suggested here or here or here.
Alternately, use your favorite alignment-based technique, if you have a good idea of what your contaminant is. In principle PB reads will be better filtered by this method than Illumina as local alignment will be able to figure out the errors quite well with longer queries.
If you go the route of trying to do alignment to some set of contaminant references, the obvious answer would be to use minimap2 as your aligner as this is exactly the task that it is designed for. You can then use samtools to filter out anything that maps well.
Here is a related question/answer that also addresses the contaminant issue and suggests using Blobology, which I'm not familiar with but may be helpful (granted it also seems to be Illumina-centric). Genome-size estimation using k-mers on PB reads appears to be pretty accurate too, so I think that it's likely ok to work with raw PB k-mers.


Setting the legChrom argument to NA will omit plotting a legend.

Found in the manual for chromPlot on Bioconductor for release 3.7.


taxon exclude list for searching local blast database using blastn

Hurray, the solution seems already to be implemented in the new BLAST database format (v5), which is usable with blast+ versions >= 2.8.0.

It is possible to specify an exclusion as well as inclusion list of taxon ids, which, however, have to be on species level.

Species level taxon ids can be obtained from higher level taxon ids using the script get_species_taxids.sh, e.g.: get_species_taxids.sh -n Enterobacterales

See docs here: https://ftp.ncbi.nlm.nih.gov/blast/db/v5/blastdbv5.pdf

Edit:
For creating taxon ID lists locally from the NCBI taxonomy files, the program subtree can be used. See example therein for compiling a list of taxon IDs from environmental or unclassified taxa.


Global or local alignment on sequences for which it is assumed that they have common ancestry

My initial reaction is that 21 amino acids is too little information to demonstrate these proteins are homologues (which is the question behind your question). A homologous site means it shares a common ancestry and is distinguishable from random chance or convergent evolution. A common ancestor can be demonstrated using phylogeny (trees).

If a regular alignment program cannot align these sequences, such as clustalo or muscle, the assumption they are homologues cannot be easily tested using phylogeny. That would then become complicated at which point you would need to ask whether it is worth pursing the question.

My recommendation is to perform a simple protein-protein blast and assess the e-scores to get some feel about you data. Your question is achievable, but you would solve this in context to additional sequence data. Whether that sequence data is a denser sample or else contiguous sequence data needs further exploration.


Raw vs Filtered in the output of cellranger count

I'm unsure whether this is the answer you are looking for, but when looking into 10X cellranger documentation for the Matrices Output:

  Unfiltered gene-barcode matrices:
  Contains every barcode from fixed list of known-good barcode sequences. This includes background and non-cellular barcodes.

  Filtered gene-barcode matrices:
  Contains only detected cellular barcodes.

Source


bedtools getfasta can do this. See documentation.


Combine structures in the same session - PyMOL

Open the PyMOL session.
Use the load command to load your PDB file(s): load file.pdb.
Use the align or super commands to overlay structures, e.g. align prot1, prot2.


Problem: "Pair-end" reads scRNA seq data (Drop-seq)

I received the following message from the author of another Drop-seq paper:

  After much conversation with the folks at the SRA, it appears the SRA
  compression algorithm automatically deletes certain subsets of bam
  tags and does so in cryptic way (in the sense of not alerting the
  uploader). Obviously this is going to be an issue for many single cell
  datasets encoded this way, and I hope the SRA will eventually change
  their system to accommodate.

If you look at the "Download" tab of the same page, there now exists a link to download the originally uploaded BAM file.

Alternately, you can use the SRA AWS tutorial to download the files on that platform.


Transferring genomic features on new coordinates

I think the standard way of doing this is to make a chain file, then use it to do a liftOver of the annotations:

mkdir psl
for i in ../ci3/rm/masked/*.masked; do blat ../ci2.2bit $i -tileSize=12 -fastMap -minIdentity=98 psl/`basename $i .fa.masked`.psl -noHead -minScore=100; done

Translate psl files to chains in the directory chain:

mkdir chain
for i in psl/*.psl; do axtChain -linearGap=medium -psl $i ../ci2.2bit ../ci3/ci3.2bit chain/`basename $i .psl`.chain; done

Merge short chains into longer ones into the directory chainMerge:

mkdir chainMerge
chainMergeSort chain/*.chain | chainSplit chainMerge stdin -lump=50

concat and sort the chains:

cat chainMerge/*.chain > all.chain
chainSort all.chain all.sorted.chain

Need info about chromosome sizes for netting:

twoBitInfo ../ci3/ci3.2bit ci3.chromInfo
twoBitInfo ../ci2.2bit ci2.chromInfo

Netting: identify alignable regions from chains:

mkdir net
chainNet all.sorted.chain ci2.chromInfo ci3.chromInfo net/all.net /dev/null

Finally, select the right alignable regions using the nets, creating a "liftOver" file:

netChainSubset net/all.net all.chain ci2ToCi3.liftOver

Run liftOver:

CrossMapy.py bed ci2ToCi3.liftOver test.hg18.bed


What does that mean if there is no bifurcations in my time series single cell data?

The most obvious interpretation is that the four cell lines are related in a linear pathway rather than a branching one.


Why I am getting different gene lengths for the same genes in different samples with rsem-calculate-expression?

If there are differences in the expression of individual transcripts between samples then both the length and effective length will differ between samples. It is commonly the case that the ratio of longer and short transcript expression varies a bit between samples, which gives rise to this.


How to separate reads, one for each reference haplotype.

Map all the reads against the reference genome
Use GATK to call SNPs
Use GATK tool [FastaAlternateReferenceMaker][1] to build the alternative reference
You now have the two "haploid references" (quotes needed because you have to take into account that they will not necessarily reflect the exact haplotypic structure, it could be a mix of the two haplotypes)
Align the reads against the two references, and attribute each read to the reference to which they map better. My colleagues developed for a different purpose (detection of Allelic Imbalance) a pipeline that includes the step for determining to which of two references a read belongs. The tool is called [BayesASE][2]

In general you may study all the work that has been done in the field of allelic imbalance and allele specific expression, where this problem has been encountered and solved.
[1]: https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker
[2]: https://github.com/McIntyre-Lab/BayesASE


After KEGG and GO analysis, how to make tables+phylogenetic trees

I'm going to answer 1+2 because 3 is hard to understand and maybe trivial?

Each of 1+2 seems to be have these sub-tasks:

identify the contigs in your assembly that have the genes of interest.
pull out the gene sequences from the contigs (you should have intermediate tblastx output somewhere with coordinates, if you don't you can just re-blast those contigs).
possibly translate those in silico to protein sequences.
tblastx (gene) or blastp (protein) those gene/protein sequences against whatever db is relevant.
pull sequences for the hits above some threshold for each gene from the db.
perform a multiple sequence alignment for each gene against the hits.
plug the MSA into whatever phylogenetic tool you like. I used phyml a lot back in the day, don't know what the kids like now.
make a visualization of the tree, probably?


How to annotate optimally a fungal genome without RNA-seq evidence?

The Maker documentation does include some instructions for training ab initio gene predictors, but it assumes an abundant EST database is available. (Assumptions about what kinds of sequences will be available for a draft genome assembly have changed drastically in the last decade.)
It may be worth exploring an iterative approach in any case. You can do a first-pass annotation with Maker using the evidence you discussed, along with ab initio predictions from Augustus (and I'd recommend SNAP as well). Each gene model will be scored by Maker with the annotation edit distance (AED), representing how well the gene model agrees with the evidence (0.0 is best). From this first pass, if you can identify a few hundred reliable gene models with good AED scores (or maybe even as many as one or two thousand), that should be plenty of training data. You could use these gene models to train a new model for Augustus, using Maker's instructions or Augustus' own instruction. (In my experience, it was quite a bit of work in either case.)


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

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

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


How to convert a .fastq chip-seq file into a wig/bigwig track and how to extract chip-seq peaks?

Before you proceed, it is usual to do quality control on the .fastq data.  The FASTQC program is a good place to start.

Once you have checked the quality, trimming is the typical next step to remove low-quality data.  There are options such as Trimomatic, trimGalore for this.

Once you are sure you have clean data, your first step is alignment - align the sequences in the .fastq file to the reference sequence in the .fasta file.  bowtie2 is one option for this process.  This gets you to the .bam stage.

These steps will take a little time, some research, and may be somewhat frustrating, but they are the foundation for what you want to do. One option for somewhat automating these steps is to use the Galaxy platform (usegalaxy.org), otherwise these are mostly command-line tools.

The bam to wig and wig to bigwig steps are relatively simple conversions, and Biostars has answers for these.


Visually summarizing phylogenetic networks

Firstly you treat this as two separate independent calculations.

The traditional approach is simply to produce a consensus phylogeny using your Bayesian method and assess the equivocal bootstraps/Bayesian probability.
Feed your data set (not the trees) into NeighborNet within SplitsTrees. There are statistics that classify the level of departure from a bifurcating phylogeny.

Examine both outputs for orthogonal support for incongruence.

Note Phylogenetic network statistics are a relatively soft area of phylogeny (they are not used as hard evidence).

Networks from Bayesian trees

The final approach (which is what you are asking about), is to feed your Bayesian trees into Dendroscope, this is a GUI that represents the uncertainty between different bifurcating phylogenies as networks and you can use it to understand the incongruence between different trees within your Bayesian MCMC, because the phylogenetic networks directly corresponds to the level of incongruence.

Dendroscope is a nice Java API and is very easy to use.


Strange p-value histogram for differential gene expression analysis

Try the protocol here for generating a different (smaller) dispersion estimate

https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues


I want to display color key on my heatmap


What is a good RNA seq normalization method that allows for across sample comparisons and between transcripts

Unfortunately, there doesn't seem to be any such metric. We use counts two-way normalized across genes and samples but even that is approximate as there is no way to ensure that all genes being counted are definitely expressed across all samples, and that all samples are sequenced to the exact same depth and that there are no batch effects and that biological differences are reflected with full fidelity.


How can I convert codon coordinates to genomic position?

HFE|LOC108783645 is https://www.ncbi.nlm.nih.gov/gene/108783645. You can use biomaRt to go from the LOCxxxx ID to location. Searching by HFE will give you a different location though, so you might have to account for such eventualities. It is a long and often manually intensive process.
You should be able to use entrez-direct/eutils to go from rsID to position
Unless there's some sort of natural language processing API, codon 109 of exon 3 of beta globin gene cannot be parsed by existing APIs (Unless it's always codon X of exon Y of Z, in which case you could use grep too get the components. Even so, going from exon + codon + codon base change to exact base location will require a lot of wrangling.


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


How to Installing & Using RDP Classifier through bioconda?

The benefit of using the conda package is that you no longer need to know or care exactly where the jar files are (that many java-based packages require this is TERRIBLE design). The example java -jar ... command in your post becomes:

Clustering derep -m '#=GC_RF' -o derep.fa all_seqs.ids all_seqs.samples alignment/*.fasta

in the bioconda version. In bioconda packages, there are wrapper (usually python) scripts meant to provide somewhat friendlier interfaces to java programs that otherwise lack executables. In this case, RDP Classfier has the following executables in bioconda:

AbundanceStats
AlignmentTools
Clustering
FrameBot
KmerFilter
ProbeMatch
ReadSeq
SeqFilters
SequenceMatch
hmmgs
classifier

If you really want to know, the .jar file is under $PREFIX/share/rdptools-$PKG_VERSION-$PKG_BUILDNUM/, where $PREFIX is where ever the environment is and the other two variables indicate the package version and build number (at the moment, rdptools-2.0.2-1 would be the likely folder name).


This indeed looks like some sort of binary file. I have seen bunch of those when people renamed .fastq.gz to .fastq without actually uncompressing them. However, as you say, this is not your case, as file pair1.fastq would identify the file as gzip compressed data.
Looking at possible return values of file comments here, data basically mean the command could not figure out anything about it beside being some sort of binary data file, but I suspect it probably somehow broken file. Honestly, at this point I think you should check the source of the data. If it's from a public database, re-download the data. If it's your own sequencing, contact the sequencing provider and check the data with them. I don't see any way to rescue this file (which does not necerasrily mean there is nothing you can do).


Pymol, in either its commercial or open source edition, will be sufficient to make figures such as your examples. Pymol recognises "polar contacts" for the intermolecular interactions and has labels that can be used for residue codes.
It is also scriptable and is widely used.


Obtaining data table headers from GEO using GEOquery

The short answer is that if you are seeing information in a GEO sample, then it is the sample that you need to download and access using the appropriate method. In this case, that method is Columns():

library(GEOquery)
gsm <- getGEO("GSM2644971")
Columns(gsm)

          Column      Description
1         ID_REF
2         VALUE       Normalized (provided the normalization method) Average Beta
3         Detection   Pval

The longer answer: you will really benefit from sitting down with the GEOquery documentation and trying to understand the objects that it uses and their methods. It is not something you can work out by trial and error for every case-by-case basis. The tutorial is a really good starting point (you will find Columns() in there, for example).


mBED has time complexity $O(n \lg n)$, claimed by Clustal Omega paper, why?

Turns out that I was dumb: I misread.
The Clustal Omega paper (in the supplementary materials) chose $t = \log N$, which I think originates from this claim from this part of the Linial et al. paper:

In random polynomial time $X$ may be embedded in $l_{2}^{O(\log n)}$ with distortion $\leq(1+\epsilon) \cdot c_{2}(X)$ (for any $\epsilon>0)$

Several questions still remain though:

The newest implementation of Clustal Omega (that I downloaded) has source code that instead chooses $t = (\log N)^2$, why changed?
The original mBED paper chooses $t = (\log N)^2$, but it was not justified on how they chose this number instead of $\log N$


It's case specific
Docking small compounds well is problematic and you need to consider whether:

if a single conformer is okay
if rigid sidechains are okay
if implicit water is okay
if a non-polarisable model is okay
if a non-covalent interaction is okay
if you are okay with a score that is a rough ΔΔGibbs approximation of the interaction (see this SO answer) or do you want something more like a K_d

Do note that for any experiment you do to have any validity you need to have controls. Do you know what can and cannot bind already? It is imperative that you score these too.
Free
I am assuming that you want to use a free program. Virtual compound screens (VCS) are a big business in pharma. Therefore, a lot of really good software is not free, such as ICM-Dock, which outscores AutoDock discuss below in most tests.
AutoDock: simple cases
AutoDock Vina is easy to use and has a PyMOL plugin (I believe), but the peptide, in both its sidechains and its backbone, is rigid, which is not good. The ligand rotamers are generated for AutoDock at runtime, but there is a hardcoded low limit, it technically isn't a rigid ligand, but nearly.
It is the best-scoring free rigid peptide docking program.
There are lots of tutorials on using Autodock, such as the official one or on YouTube. As shown by this Bioinformatics Stack Exchange question Autodock 4 (different score fxn) has several steps, each important.
Rosetta: very complicated cases
I prefer the harder to use Rosetta Ligand_dock which uses conformers (if the ligand is parameterised correctly, say they were generate via Open Babel conformers) and allows sidechains to repack.
Water
Both tools mentioned use implicit waters, which works badly for highly hydrophobic ligands.
That is the waters are simulated as a homogeneous field that affects the ligand and protein. But some solutions are possible, in the case of AutoDock WaterDock variant exists —although this is not polarisable water.
For more advanced cases there are MD methods such as dynamic undocking, which are completely non-trivial. All these methods mentioned use classical force-field calculations (Amber, Charmm, Talaris etc. are different models), but in MD noise (think of it as environmental heat) is present in the calculations, so the ligand comes unbound (and is actually pulled out). You can also get better results by refining your bound ligand with QM-MD (e.g. gauss of gaussian).
Choices...
This is a lot of options. So a first question is are your ligands man-made drugs/fragments (i.e. obey Lipinski's rule of 5) or metabolites? Generally the former are rigid and very hydrophobic (greasy) (and require better force fields), while the latter are flexible and hydrophilic.
Minor considerations
Also to consider is your entranceway to your protein as the Michaelis constant is dictated in some enzymes by the substrate going in, not by its binding (e.g. P450s).
If the ligand have a shared backbone that is fixed there are some tricks you can do to cut your compute time and increase the reliability.
If the ligands are bound not all methods are possible. Rosetta ligand_dock can but it is a bit fiddly to parameterize.
RDKit skills
One thing also to consider is your programming skills. Can you use RDKit? In most cases downloading the sdf off PubChem and running OBabel or similar to generate conformers is not possible and you have be familiar with a python package called RDKit, which is really powerful but non-trivial to use.
Additionally, docking scores can be strongly augmented with ligand-only values (such as logP, TPSA, molar refractivity, molecular weight, QED, Bertz topology etc. etc) that are generated from RDKit.
NB. The scoring function in Autodock Vina (not 4 though, as that is Amber FF) partially takes into account.
Ligand_dock details
Ligand_dock is an older application of the Rosetta suite and the current philosophy is to use a specific Rosetta script, however, it works fine and is actually quicker to implement. Here are some links:

manual of ligand_dock
tutorial for ligand parameterisation
PubChem for your SDF files
online tools for small molecule manipulation
example of ligand creation for docking with RDKit


10x Genomics Chromium single-cell RNA-seq data analysis options?

Data preparation

Cell Ranger uses the Illumina sequencing output (.bcl) files

Make fastq files:
cellranger mkfastq ==> .fastq
Prepare count matrix: cellranger count ==> matrix.mtx, web_summary.html, cloupe.cloupe
Optional: combine multiple matrix.mtx files (libraries): cellranger aggr

Data analysis

Loupe Cell Browser visualization of cloupe.cloupe files

Count table matrix.mtx analysis options:

Python
R Cell Ranger R Kit: cellrangerRkit::load_cellranger_matrix() ==> ExpressionSet
R Scater: scater::read10XResults() ==> SCESet object
R Seurat: Seurat::Read10X() ==> Seurat object


Smith-Waterman is guaranteed to give you the best alignment given your choice of substitution matrix.  For a few such tiny sequences, it's what I would use.  Many implementations exist.


STRINGdb: reverse map STRING identifiers to gene names?

This answer does not rely on STRINGdb package.
STRING provides supplementary identifier mapping files here: https://string-db.org/mapping_files/.
It looks like you are interested in human STRING name one - https://string-db.org/mapping_files/STRING_display_names/human.name_2_string.tsv.gz. You could read it as a separate lookup.


Phylogenetics software that can represent collapsed nodes as triangles

For Mac Figtree here is highly recommended, it is also available for Windows and Linux here. The technique is called "collapsing clades" and is the only package I know that does it. I would recommend using it, because I would assume it is averaging the branch length of a given clade and the 'base' of the triangle is proportional to the number of taxa. Thus it recognised as being very reliable, albeit the precise algorhithm it uses ... I don't know :-O

With FigTree you get two different "triangle options" in the menu, which give different shaped triangles for collapsing the number of lineages within the select clade you chose. The big one, is the algorhithm appropriately described above. It has many options for producing publication quality diagrams, albeit I then import the output into a graphics program for further editing (there are shortcomings IMO in colouring and fonts. Then again for free is better than good.

If you are outputing a tree within a script, ETE3 (Python) will do a "collapsing of clade" option, but it isn't pretty triangle stuff. ETE3 will produce a publication quality tree diagram from this and reduce to 'low resolution', which is very useful in large analyses of loads of trees.

In summary, I don't think FigTree has any competition in the context you are seeking.

Moreover, you can configure your bin, or alias file so you can automatically open a treefile just typing,

figtree treefile.tre

Also you can reduce you font size to next point 1 or 2 (okay you can do that in anything), if you are attempting optical analysis on loads of taxa. Its also really stable. I just find it really useful, like a "jazz standard" its a "phylogeny standard".


Papers

The space-efficient BWT construction algorithm is in bwt_gen.c. One of the first things you see there is the copyright line:

  Copyright (C) 2004, Wong Chi Kwong.

With that name, you can find the relevant papers:

Lam et al: Compressed indexing and local alignment of DNA (Bioinformatics, 2008). This paper describes the BWT-SW alignment algorithm.
Hon et al: Practical Aspects of Compressed Suffix Arrays and FM-index in Searching DNA Sequences (ALENEX 2004). This was the first paper to demonstrate that it's feasible to build the FM-index for a human genome.
Hon et al: A Space and Time Efficient Algorithm for Constructing Compressed Suffix Arrays (Algorithmica, 2007). This paper contains the theoretical description of the construction algorithm.

The third paper only mentions BWT briefly in the last section. It deals with compressed suffix arrays, which are space-efficient representations of the inverse function of LF-mapping. But if you have a representation of (the inverse of) the LF-mapping, you have a representation of the BWT. Hence the algorithm can be adapted for BWT construction.

Algorithm

When building BWT, you usually append a unique character $ to the end of the text, either implicitly or explicitly. This way you can reduce sorting rotations to sorting suffixes. You sort the suffixes in lexicographic order, and then you output the preceding character for each suffix in that order to form the BWT.

Assume that you already have the BWT for string X, and you want to transform it into the BWT of string cX, where c is a single character. You can do that in three steps (this assumes 0-based array indexing):

Find the position i of the full suffix X, and replace the $ in the BWT with the inserted character c (set BWT[i] = c).
Determine the number of suffixes j that are lexicographically smaller than cX with j = LF(i, c) = C[c] + BWT.rank(i, c). Here C[c] is the number of suffixes that start with a character smaller than c, and BWT.rank(i, c) is the number of occurrences of c in the prefix BWT[0..i-1].
Insert $ between BWT[j-1] and BWT[j].

By iterating this, you can build the BWT for any string slowly but space-efficiently. The algorithm of Hon et al. is based on the same idea, but they make it faster by inserting multiple characters at the same time. Given the BWT of string X and another string Y, they produce the BWT of the concatenation YX. The details are quite complicated, but you can find them in the third paper and in the source code.


MiXCR: only create a single export file for all clonotypes

This is from the MiXCR documentation:

Under the hood analyze amplicon is equivalent to the execution of the following
MiXCR actions options:
> mixcr align -s hs -p kAligner2 input_R1.fastq input_R2.fastq alignments.vdjca
... Building alignments
> mixcr assemble alignments.vdjca clones.clns
... Assembling clones
> mixcr exportClones --chains IGH clones.clns clones.txt
... Exporting clones to tab-delimited file

So when using analyze amplicon, you are using mixcr exportClones under the hood and not exportAlignments, hence the difference in output.

I believe what you are experiencing results from a bug, please see this issue. If this is the case, all these files should be the same according to the bug report. Fortunately the developers have already addressed the bug in version 3.0.3.


polishing assembled genome to QV50 value

15 times is too many times. If you were going to get quiver to get you good QVs it would have happened after 2-3 rounds probably.

Some questions:

Why are you using quiver? Are your reads from RSII chemistry? Arrow is recommended for Sequel reads.
What is your dataset coverage? If you do not have very high coverage (I think >30X is usually a minimum for a eukaryotic genome) it may be that you are just not going to get exactly the values you want because assembling low coverage data is hard.
Have you tried other polishing tools? They might be able to get you better value. I am not an expert but some people like to use Racon or Apollo.
While I'm not an expert, looking into how the QVs are calculated can be helpful. I looked at the quiver docs and it looks like if the reads map ambiguously then QV is low. So if you for example have lots of haplotigs, then mapping might be ambiguous for a lot of reads, and therefore this will depress your QV values. This is therefore a problem that can come from the starting assembly and not just the polishing. You may want to go back and take another look at your starting assembly.


You can do this easily with bioawk, which is a version of awk with added features facilitating bioinformatics:

bioawk -c fastx '{print $name"\t0\t"length($seq)}' test.fa

-c fastx tells the program that the data should be parsed as fasta or fastq format. This makes the $name and $seq variables available in the awk commands.


How to extract hypervariable region of mitochondrial dna from its fasta file?

You can use Mitomap to get the coordinates of the hypervariable regions in the rCRS mitochondrial sequence (the standard MT sequence in the field). I checked the Genome Loci section of Mitomap and saw that it lists three:

Map Locus   Starting    Ending  Shorthand   Description
MT-HV1      16024       16383   CR:HVS1/HV1 Hypervariable segment 1
                                            [classic:16024-16365]
MT-HV2      57          372     CR:HVS2     Hypervariable segment 2
                                            [classic:73-340 ]
MT-HV3      438         574     CR:HVS3     Hypervariable segment 3

That gives you the coordinates on the rCRS sequence. What you can do now, is align all of your sequences and the rCRS together and see if you can infer where the hypervariable regions fall in your target sequences. Since, presumably, the rest of the sequence will be more conserved since it isn't hypervariable, you should be able to identify the syntenic blocks and locate the hypervariable regions in your sequences. I would expect to see a very high level of sequence identity in the non-hypervariable regions across all of your sequences (since they're all human mitochondria), so you should be able to find the hypervariable ones as the regions of low sequence similarity which also coincide with the known coordinates of the rCRS.


What unit (TPM, FPKM/RPKM, or other) to use when working across samples

Neither of RPKM, FPKM or TPM is a good choice. All these methods are similar, they only perform the operations in a slightly different order. Eventually they correct for gene length and sequencing depth. Correction for sequencing depth is necessary but not sufficient in many cases. Correction for length is typically not done for a differential analysis as it reduces the counts of many genes and therefore throws away statistical power. I will not describe this in detail since the StatQuest video series you link includes videos about DESeq2 and edgeR normalization procedures which extensively cover the normalization procedures and why these approaches are superior to RPKM/FPKM/TPM. In short: You have to correct for differences in library composition to compare between samples. Please watch these videos and read the manuals of the established RNA-seq analysis packages such as but not exclusively limma, edgeR and DESeq2.


What motif finding software is available for multiple sequences ~10Kb?

The MEME Suite web site contains a collection of tools for motif analysis (I'm one of the maintainers). It contains two de novo motif discovery tools: MEME and DREME. Public web applications are provided, but you can also download and build command line tools for a local installation.

For your first goal you could use MEME and select the "Any number of repetitions model" (ANR). For your second goal, you'd use MEME with the "Zero or One Occurrences Per Sequence" (ZOOPS) model. For your third goal you could use FIMO (Find Individual Motif Occurrences), and one or more of the motif databases provided on the software and database download page.

It sounds like your sequence data is about 30Mb. The MEME web application is limited to 60kb of sequence data, so you'd have to install a local copy of the MEME Suite. MEME would take a long time to analyze a 30Mb sequence database unless you have MPI configured, and lots of cores available. You might want to consider analyzing a randomly selected subset of your sequences. The running time of MEME grows as the cube of the number of sequences.

For short motifs, you may want to use DREME rather than MEME. DREME is better than MEME at identifying short motifs, but is limited to motifs <= 8 positions wide.


Is there a command line tool to split a SAM/BAM file by CB (cell barcode) tag?

I wrote a python script based on pysam that is free for anyone to use:

"""
MIT License
Copyright (c) 2020 Warren W. Kretzschmar
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

from pathlib import Path

import pysam

class BamWriter:
    def __init__(self, alignment, barcodes, prefix):
        self.alignment = alignment
        self.prefix = prefix
        self.barcodes = set(barcodes)
        self._out_files = {}

    def write_record_to_barcode(self, rec, barcode):
        if barcode not in self.barcodes:
            return
        if barcode not in self._out_files:
            self._open_file_for_barcode(barcode)
        self._out_files[barcode].write(rec)

    def _open_file_for_barcode(self, barcode):
        self._out_files[barcode] = pysam.AlignmentFile(
            f"{self.prefix}_{barcode}.bam", "wb", template=self.alignment
        )

def main(input_bam, barcodes_file, output_prefix, contigs):
    """Split a 10x barcoded sequencing file into barcode-specific BAMs

    input:
    barcodes_file: a file containing barcodes or a single barcode
    contigs: '.' for all contigs, 'chr1' for the contig 'chr1',
    or '1-5' for chromosomes 1, 2, 3, 4, and 5
    """
    alignment = pysam.AlignmentFile(input_bam)
    if Path(barcodes_file).is_file():
        with open(barcodes_file, "r") as fh:
            barcodes = [l.rstrip() for l in fh.readlines()]
    else:
        barcodes = [barcodes_file]
        print(f"Extracting single barcode: {barcodes}")
    writer = BamWriter(alignment=alignment, barcodes=barcodes, prefix=output_prefix)
    if contigs == ".":
        print("Extracting reads from all contigs")
        recs = alignment.fetch()
    else:
        if "-" in contigs:
            start, end = contigs.split("-")
            print(f"Extracting reads from contigs {start} to {end}")
            recs = (alignment.fetch(str(contig)) for contig in range(start, end + 1))
        elif "," in contigs:
            contigs = contigs.split(",")
            print(f"Extracting reads from contigs {contigs}")
            recs = (alignment.fetch(str(contig)) for contig in contigs)
        else:
            print("Extracting reads for one contig: {contigs}")
            recs = (alignment.fetch(c) for c in [contigs])

    for rec in recs:
        try:
            barcode = rec.get_tag("CB")
            writer.write_record_to_barcode(rec=rec, barcode=barcode)
        except KeyError:
            pass

if __name__ == "__main__":
    import sys

    main(
        input_bam=sys.argv[1],
        barcodes_file=sys.argv[2],
        output_prefix=sys.argv[3],
        contigs=sys.argv[4],
    )
```


Mini-bulk RNA-seq and scRNA-seq analysis differences?

Single-cell tools were developed to address some specific factors unique to scRNA-seq, such as:

a lot of dropouts, so hundreds or thousands of detected genes, not tens of thousands (for mammalian genomes)
hundreds or thousands of replicates
not necessarily well-defined groups

Thus, I would say your experiment is a lot closer to bulk RNA-seq than scRNA-seq.


How to check whether all BAM read contain defined read groups?

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

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

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

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

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

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


How to superimpose barplots in karyotype graphic using the chromPlot Bioconductor package?

chromPlot(gaps = cel_gap, figCols = 2, annot1 = annot1, annot2 = annot3, annot3 = annot2,
          chrSide=c(1, 1, -1, 1, -1, 1, -1, 1))

The chrSide argument decides on which side the histogram goes, so it can be used to solve my problem. The bars are not stacked but the result is the same: more accurate visualisation.

Thanks to Llopis for the guidance.


How do I perform a pathway related grouping of genes?

The DAVID tool is always interesting, and easy to try: https://david.ncifcrf.gov/gene2gene.jsp

You submit a list of genes, it uses several different annotation databases to annotate the gene list, and then it clusters the results and tests for statistical enrichment. You can potentially learn about functional enrichment in your gene list from several different independent annotations.