How to convert a phylogeny to a dendrogram in R

For your example to work you should replace the rtree function with the rcoal function, see example below.

library(ape)
library(dendextend)
Tree <- rtree(10, rooted=F)
is.ultrametric(Tree)
[1] FALSE
is.binary.tree(Tree)
[1] TRUE

Tree_rcoal <- rcoal(10, rooted=F)
is.ultrametric(Tree_rcoal)
[1] TRUE
is.binary.tree(Tree_rcoal)
[1] TRUE

Tree_rcoal2 <- rcoal(10, rooted=F)

tanglegram(Tree_rcoal, Tree_rcoal2)

To make use of a file in newick format, you can try to convert with the chronos function.

s <- "owls(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3,Tyto_alba:13.5);"
cat(s, file = "ex.tre", sep = "\n")
tree.owls <- read.tree("ex.tre")
is.ultrametric(tree.owls)
[1] FALSE
dendrogram <- chronos(tree.owls)

Setting initial dates...
Fitting in progress... get a first set of estimates
         Penalised log-lik = -24.04992
Optimising rates... dates... -24.04992
Optimising rates... dates... -24.04908

Done.
is.ultrametric(dendrogram)
[1] TRUE


Reading fasta sequence regions with faidx

If you can't read in the entire sequence for a given reference genome:
(1) Use a script to generate a BED file of masked regions of low complexity (i.e. where you find Ns in your FASTA file). (2) Use bedops --complement with this mask file and the output of UCSC fetchChromSizes to generate a BED file containing regions that are unmasked. (3) Repeat for other reference genomes, as needed.
You can now loop over the unmasked regions in this file, making calls to samtools faidx or faidx_fetch_seq(). These calls will return unmasked sequence.
For performance or convenience, you could read chunks of sequence into a std::deque or std::vector buffer and then read through that however many hundreds of bases, refreshing the buffer, as needed.


I'm one of the authors of HiGlass. Sorry you're getting that error. I just tried it out in a fresh environment and it worked for me. Just one thing is missing in the docs which is that you need to do import higlass.

Here's what I did to set up my environment:

conda create --name test-bi python=3.6
source activate test-bi

The rest is from the docs:

pip install jupyter higlass-python
jupyter nbextension install --py --sys-prefix --symlink higlass
jupyter nbextension enable --py --sys-prefix higlass
jupyter notebook

It looks like either you or somebody else added a bug report here: https://github.com/higlass/higlass-python/issues/29

I'll respond and continue the conversation there.

Also, there's a repository full of examples from the HiGlass SciPy presentation here: https://github.com/higlass/scipy19


How to convert a binary matrix of gene presence or absence into a fasta sequence

It's not a fasta file, but:

> m
       sample1 sample2 sample3 sample4
fliI         1       1       1       1
patB_1       1       1       1       1
pgpA         1       1       1       1
osmB         1       1       1       1
cspA         1       0       1       1
> # Collapse to labeled strings
> blah = apply(m, 2, function(x) paste(x, collapse=''))
> blah
sample1 sample2 sample3 sample4
"11111" "11110" "11111" "11111"
> # Write that to a file with the appropriate format
> cat(paste(mapply(function(x, y) sprintf(">%s\n%s\n", x, y), names(blah), blah), collapse=""), file="some file.fa")

Make sure to change some file.fa. I have no idea why you would want to do all of this, but this will give you the output you want:

$ cat blah.txt
>sample1
11111
>sample2
11110
>sample3
11111
>sample4
11111


GenomicRanges: get nearest neighbor distances for random genes using a for loop

Is it because rand.gr has not been initialised, so when you ask to set rand.gr[i] it can't because there isn't one yet?

If you set i to a value (e.g. 1), and run the lines inside your loop, which of the lines produces the error?

If you don't need to retain all the random genomic ranges at the end, you could just leave out some of the [i]s, something like this:

for(i in 1:ncol(rand.df)){
  rand.gr <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists <- distanceToNearest(x = rand.gr)
  output <- print(mcols(dists[,1]))
}


How to extract all pair end fastq.gz files from multiple subdirectory to a single folder

Go to the common parent folder under which all FASTQ.GZ files reside. They could be in sub-folders under the folder, but all of them should have the folder you're in as a common parent folder. From such a folder, you should be able to run a find

find . -name "*.fastq.gz"

to get your FASTQ.GZ files. Now, there are surely better ways of doing this but here's one straightforward way.

Let's say you need to extract the FASTQs to ~/target/dir/, you'd do something like:

for f in $(find . -name "*.fastq.gz")
do
    tgt_file="~/target/dir/${f%.*}" #removed .gz from the name
    gzip -dc $f >$tgt_file #pigz -dc would be better if you have pigz
done

EDIT: OP has clarified that they wish to copy the files as such, not decompress them.

There are two ways to achieve this goal - you can either make physical copies of the files (thereby doubling the storage space you're currently using for them), or you can create symbolic links to the files (which are like shortcuts pointing to the actual files). I'll show you both ways here:

Let's assume the files are all under /fastq/files

find /fastq/files/ -name "*.fastq.gz" -exec cp {} ~/target/dir/ \; #This creates physical copies

cd ~/target/dir
find /fastq/files/ -name "*.fastq.gz" | xargs -I v_f ln -s v_f #These lines create symbolic links that point to the original files. Do not delete/move/rename the original files or these will become invalid.


How to calculate the memory usage of storing kmers in RAM

To follow up on Devon Ryan's answer, I thought it would be a little fun to write a Python script that demonstrates using a bit array to maintain a presence/absence table.

Note: I wrote a C++ port that includes a custom bitset implementation that can be sized at runtime. This and the Python script are available on Github: https://github.com/alexpreynolds/kmer-boolean

#!/usr/bin/env python

import sys
import os
import bitarray

# read FASTA
def read_sequences():
    global seqs
    seqs = []
    seq = ""
    for line in sys.stdin:
        if line.startswith('>'):
            if len(seq) > 0:
                seqs.append(seq)
                seq = ""
        else:
            seq += line.strip()
    seqs.append(seq)

# build and initialize bit array
def initialize_bitarray():
    global ba
    ba = bitarray.bitarray(4**k)
    ba.setall(False)
    sys.stderr.write("Memory usage of bitarray.bitarray instance is [%ld] bytes\n" % (ba.buffer_info()[1]))

# process sequences
def process_sequences():
    global observed_kmers
    observed_kmers = {}
    for seq in seqs:
        for i in range(0, len(seq)):
            kmer = seq[i:i+k]
            if len(kmer) == k:
                observed_kmers[kmer] = None
                idx = 0
                for j in range(k-1, -1, -1):
                    idx += 4**(k-j-1) * bm[kmer[j]]
                ba[idx] = True

def test_bitarray():
    test_idx = 0
    for j in range(k-1, -1, -1):
        test_idx += 4**(k-j-1) * bm[test_kmer[j]]
    test_result = ba[test_idx]
    if test_result:
        sys.stdout.write("%s found\n" % (test_kmer))
        sys.exit(os.EX_OK)
    else:
        sys.stdout.write("%s not found\n" % (test_kmer))
        sys.exit(os.EX_DATAERR)

def main():
    global k
    k = int(sys.argv[1])

    global bm
    bm = { 'A' : 0, 'C' : 1, 'T' : 2, 'G' : 3 }

    read_sequences()
    initialize_bitarray()
    process_sequences()

    try:
        global test_kmer
        test_kmer = sys.argv[2]
        if len(test_kmer) == k:
            test_bitarray()
        else:
            raise ValueError("test kmer (%s) should be of length k (%d)" % (test_kmer, k))
    except IndexError as err:
        keys = list(observed_kmers.keys())
        for i in range(0, len(keys)):
            sys.stdout.write("%s found\n" % (keys[i]))

    sys.exit(os.EX_OK)

if __name__== "__main__":
    main()

Note that this doesn't look at canonical kmers, e.g., AG is considered a distinct 2mer from its reverse complement CT.

To use this script, you pipe in your FASTA, specify the k, and an optional kmer that you want to test for presence/absence, e.g.:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 2 CG
Memory usage of bitarray.bitarray instance is [2] bytes
CG found

Or:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 3 AAA
Memory usage of bitarray.bitarray instance is [8] bytes
AAA not found

Or if the optional test kmer is left out:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 5
Memory usage of bitarray.bitarray instance is [128] bytes
CATTC found
ATTCT found
TTCTC found
TCTCG found
CTCGG found
TCGGG found
CGGGA found
GGGAC found
TTATA found
TATAT found
TTTAT found
TTATT found
TATTA found
ATTAG found
TTAGA found
TAGAC found
AGACC found
GACCT found
ACCTC found
CCTCT found

Or for the ~67M kmers in a 13-mer set, for which a roughly 8.4MB bit array is reserved:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 13
Memory usage of bitarray.bitarray instance is [8388608] bytes
TTTATTAGACCTC found
TTATTAGACCTCT found


How to interpret track height in Chipseq?

Depends how the files were made. In the simplest case yes, the height represents the pipeup of reads from the BAM file that was used. For a direct (visual) comparison you have to normalize the files though as otherwise sequencing depth confounds the height of the peaks, like if file A has ten times more reads and the peak is 10 times higher than in B then this is probably only due to differences in total read number (=sequencing depth). I wrote down my preferred method to do that in this biostars post: https://www.biostars.org/p/413626/#414440
Does this make sense to you?


I think this is a bug in the input data.
At a wild guess what I think is happening is the input data is not regarded as paired, hence the error paired = paired. The other weird error is x must be numeric. Firstly,

how are you scoring missing data, its not via - by any chance? Missing data must be via nan type syntax.
are you sure there are no heading bugs? Like a header is getting dragged into the numeric data?

Wilcoxan is a non-parametric T test, so the wet root and dry root must match exactly the same sample. Thus it is not a general wet root versus generic dry root, but wet and dry from the same specimen. If this is the data you have the program doesn't think it is matching.
If it is truely paired, again I think the source of error is you have some missing value which is leaving the pairings unequal lengths or else a header bug.
If the data genuinely isn't paired, then I would look at the Mann-Whitney U test, it would be easier in this scenario.
I don't use R syntax (just pandas dataframes) so can't say whether the syntax in the first example is valid, but it is pretty clear regarding your input. In Python you'd generally just give it the entire dataframe and tell it to go figure.
The reality is that loads of programs do Wilcoxan, you don't need to use R if its a one-off test. Obviously if this a part of a pipeline (automation) thats a bit different.


What is a quick way to find the reverse complement in bash

Thanks to Manu Tamminen for this solution:

echo ACCTTGAAA | tr ACGTacgt TGCAtgca | rev


How do I computationally predict the binding position of a DNA binding protein on the promoters of differentially expressed genes?

You could run FIMO on your entire genome for TFs (transcription factors; DNA-binding proteins) of interest, which gives you binding sites: genomic intervals where those TFs bind.

https://bioinformatics.stackexchange.com/a/2491/776

You can then map or intersect promoter regions with those TF binding sites with BEDOPS bedmap or bedops, respectively:

https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html

The use of bedmap is recommended to assign TFs to a promoter of interest. This tool will return both the promoter region and any TFs that overlap it by the specified overlap criteria.


Error: sort order of reads in BAMs must be the same

I think cuffdiff is not quite giving you the right error...it looks like the crux of your problem is you aligned each sample to a different reference.


How can HISAT2/StringTie report decimal coverage values

I suspect this is average coverage across the transcript or exon.

Formula for average coverage:

Where:

N: number of reads mapping/aligning to your transcript/exon
L: average read length
T: Transcript/exon length


One option is to use one of the apply family functions. If you enclose your code above (which generates errors by the way) in a function, you can then "apply" this function on the rows or columns of your data with something like apply(X, MARGIN, FUN, …), the margin being 1 for rows or 2 for columns (as in your case).


Assuming your user is in the group permitted to use docker:

docker pull quay.io/biocontainers/tophat:2.1.1--py27_3
nohup docker run -v /some/path:/data quay.io/biocontainers/tophat:2.1.1--py27_3 tophat2 ... &

The only thing to remember is to mount everything you'll need into the docker container (/data in the example above). If you need to debug things then you can add the -it flag to docker run and use bash as the command.

As an aside, you might find using screen or tmux preferable to nohup ... &, since you can reconnect to sessions with them.


Finding new drug-protein targets/interactions using a known one as base

If you have information on which part of the protein are interacting, and its sequence, you can search for homologous or similar sequences in other proteins. This would give you a list of potential targets and probably a similarity score which you can also use. Then you would have something to run simulations on.

And I think that with 3D structure you should be able to run simulations in silico in this case (eventhough it is probably more common to look for the ligands to the protein, not proteins to the ligand).

In silico simulations are computationally intensive, but I´ve also seen some being run on a laptop, so you don´t necessarily need access to HPC-clusters, or supercomputer.


Welcome to the world of all possible pairwise comparisons in DESeq2 !

So the easiest way is to create another group in your data frame, first I simulate some sensible counts:

counts = counts(makeExampleDESeqDataSet(m=48))
design = expand.grid(id=1:3,
phenotype=rep(c("sick","healthy"),each=2),
type = rep(c("A","B"),2))

And you make a group that is a combination of phenotype and type:

design$group = paste(design$phenotype,design$type,sep="_")

Run DESeq2, including id:

dds <- DESeqDataSetFromMatrix(countData=counts, colData=design,
design = ~ id + group)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("group", "healthy_A", "sick_A"))

So this compares only healthy_A with sick_A.. hopefully this is what you need... The other option (normal statistical) is to fit a model with interaction, but I think this is easier to explain and interpret.

Another question you might have is why not do two separate analysis for A and B, I think it's sometimes better if you model the variance with all available data, given the limited number of replicates.


ATAC seq peak annotation : Which peak to choose for downstream analysis

You can see from your table that the different peaks correspond to different promoters for different transcripts, for example, FAM41C has ENST00000432963.1, ENST00000446136.1, ENST00000635557.1. And the fact that they are all within the TSS (see distanceToTSS) suggest these peaks are actually falling onto actually TSS.
The underlying issue is are these TSS or do will we mis classify the peaks that are inside genes or other non-coding elements as TSS peaks. Nothing to do with TF, histones etc.
I cannot see the number of reads in those peaks, but if one of them is really strong, and the rest almost non-existent it would make sense.
How can we improve this, I think it depends on the aim of your analysis. Since you are looking at TSS, I guess it makes sense to be able to link it to something that is actually expressed. You can do:

To use an updated annotation (i see that some of the transcripts come from old hg37 annotation)

you have prior knowledge which transcripts might be expressed from RNA seq data, I would use that to reduce the number of transcripts


Compare clusters from different datasets

I've been working with the multi-CCA alignment workflow in the Seurat package and am very impressed with functionality. As usual, the CCA tutorial is excellent for an alignment workflow and will get you off the ground quickly. You can even use the DESeq2 algorithm for determining differential expression from within Seurat.

When you use the function RunCCA() you can pass rescale.groups=TRUE if you wish to re-scale.

Abstracted workflow

Data preprocessing and gene selection
Define a shared correlation space with canonical correlation analysis
Remove dataset specific cells
Align correlated subspaces using dynamic time warping
Integrated analysis across datasets (clustering, trajectory building, etc ...)

Butler, A., and Satija, R. (2017). Integrated analysis of single cell transcriptomic data across conditions, technologies, and species.


What are some good sources of 3rd order gene expression data

I familiar with the inflammatory bowel disease  which is quite complicated:
It affects multiple sites and depending on when it appears it is different. So you can look at the GEO for microarrays and RNA-seq data of this disease and I'm sure you'll find datasets from the same patients from several tissues or time.

In this disease there is a multi'omics project that fits your needs. You can find (and download) data of high dimensionality: there are samples for the same time, tissue (but different locations), gene, over a period of time and for several omics (not only expression data).
The only problem is that there are some mistakes on the metadata.

However, I don't think they are usually described as 3rd order tensors and I doubt you'll find enough data to train a machine learning method.


Yes, there is a tool to do this called R2C2 by the Vollmers lab at UC Santa Cruz.

https://www.biorxiv.org/content/early/2018/06/04/338020


Yes, this can be done, particularly if your concern is with the clustering / UMAP rather than the transformation process; I've done this recently with some Rhapsody single cell data that we generated. The main trick is adding 'return.only.var.genes = FALSE', as well as calculating a variable set that is a bit larger than the desired set.
Note that UMAP is a technique to place cells, rather than genes. If you have the information on where the cells should go, then plotting gene expression on top of that via FeaturePlot should work without any reprojection needed.
Here's what I did (with experiment-specific variables excluded):
sct.data <- SCTransform(sct.data, vars.to.regress="percent.mt",
                        variable.features.n = 3000 + length(genesToExclude),
                        return.only.var.genes = FALSE);

sct.data <- RunPCA(sct.data, verbose = TRUE,
                   features=setdiff(VariableFeatures(sct.data), genesToExclude));
sct.data <- RunUMAP(sct.data);
sct.data <- FindNeighbors(sct.data);
sct.data <- FindClusters(sct.data);

All transformed genes should be available for querying in the results.


The paper doesn't seem to have been written for general readers. Unfortunately I don't understand the equations enough to be able to explain what the categorical and continuous variables are. Just in case it helps, here are the equations from the paper showing the precision/recall calculations:

$$\begin{array}{@{}rcl@{}} \text{pr}(\tau) &=& \frac{1}{m(\tau)}\sum\limits_{i=1}^{m(\tau)} \frac{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in P_{i}(\tau) \wedge f \in T_{i}\right)}{\sum_{f} \unicode{x1D7D9}\left(f \in P_{i}(\tau) \right)},\\ \text{rc}(\tau) &=& \frac{1}{n_{e}}\sum\limits_{i=1}^{n_{e}} \frac{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in {P}_{i}(\tau) \wedge f \in T_{i}\right)}{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in T_{i} \right)}, \\ F_{\max} &=& \max_{\tau} \left\{ \frac{2\cdot \text{pr}(\tau)\cdot \text{rc}(\tau)}{\text{pr}(\tau) + \text{rc}(\tau)} \right\}, \end{array} $$

  where P i (τ) denotes the set of terms that have predicted scores greater than or equal to τ for a protein sequence i, T i denotes the corresponding ground-truth set of terms for that sequence, m(τ) is the number of sequences with at least one predicted score greater than or equal to τ, $\unicode{x1D7D9}\left (\cdot \right)$ is an indicator function, and n e is the number of targets used in a particular mode of evaluation.


Get gene names corresponding to CPG probes

Here's a Python script (e.g., called map_entrez_to_ensembl.py that uses mygene.info to convert a list of Entrez IDs to Ensembl gene names:

#!/usr/bin/env python

import sys
import mygene

mg = mygene.MyGeneInfo()

genes = ['3815',
         '3816',
         '2341']

results = mg.querymany(genes, scopes=['entrezgene'], fields=['ensembl'], species='human', verbose=False)

for res in results:
    entrez_id = res['query']
    if 'ensembl' in res and 'gene' in res['ensembl']:
        sys.stdout.write('%s\t%s\n' % (entrez_id, res['ensembl']['gene']))

Result:

$ ./map_entrez_to_ensembl.py
3815    ENSG00000157404
3816    ENSG00000167748

Looks like there's an R-based interface to mygene.info, if you really need to use R.

You could use their API test page to see how data are structured and what keywords are available for different kinds of searches.

For instance, if you wanted aliases, in addition to Ensembl names:

#!/usr/bin/env python

import sys
import mygene

mg = mygene.MyGeneInfo()

genes = ['3815',
         '3816',
         '2341']

results = mg.querymany(genes, scopes=['entrezgene'], fields=['ensembl', 'alias'], species='human', verbose=False)

for res in results:
    entrez_id = res['query']
    if 'ensembl' in res and 'gene' in res['ensembl'] and 'alias' in res:
        for alias in res['alias']:
            sys.stdout.write('%s\t%s\t%s\n' % (entrez_id, alias, res['ensembl']['gene']))

Then:

$ ./map_entrez_to_alias_and_ensembl.py
3815    C-Kit   ENSG00000157404
3815    CD117   ENSG00000157404
3815    PBT     ENSG00000157404
3815    SCFR    ENSG00000157404
3816    KLKR    ENSG00000167748
3816    Klk6    ENSG00000167748
3816    hK1     ENSG00000167748


Annotatation of transposable elements / transposons / SINE sequences in human

repeatmasker is meant to annotate/mask a genome given a library of known transposable elements. From the project description, I assume you have not done any sequencing on your own. Then, instead of annotating, I think you should find already computed annotations of TEs.
The alternative to subscription-based repbase is dfam. And Just looking at their webpage they seem to have an annotation of all individual TEs in hg38. I bet my shoes there will be plenty of TE annotations in human in genome databases too.


Searching for proteins by primary and secondary structure

Your request is rather specific, so there is no tool...

It is a straightforwards request though. I'd impose the requirement that the structures be solved and no predicted secondary structure. You can do it two ways if using Python.

PDB

After downloading all PDB structures, in python with PyMOL (installable with conda, tricky without) you'd do something like

import pymol2, os
from multiprocessing import Pool

def find_motif(file):
    motif = 'AFG'
    try:
        with pymol2.PyMOL() as pymol
            pymol.cmd.load(file)
            seq = pymol.cmd.get_fastastr()
            if motif not in seq:
                return None
            d = {x: []}
            # str.find will cause issues if there is a gap also structure may not start at one.
            pymol.cmd.iterate('resn ALA+PHE+GLY', 'x.append((resi, resn, ss))', space=d)
            sequence = -1
            i = 0
            for residue in d['x']:
                if residue['resn'] = 'ALA':
                    sequence = 0 # lazy enum
                    i = int(residue['resi'])
                elif sequence != -1 :
                    if int(residue['resi']) == i+1 and residue['resn'] = 'PHE':
                        sequence = 1
                    elif int(residue['resi']) == i+2 and residue['resn'] = 'GLY':
                        return (file, i)
                    else:
                        sequence = -1
    finally:
        return (file, None)

if __name__ == '__main__':
    p = Pool(12)
    structure_files = [f for f in os.listdir('pdb') if os.path.splitext(f)[1] == '.pdb']
    verdicts = p.map(finder, structure_files)
    print([(f,v) for f,v in verdicts if v is not None])

The end result will be a list of PDB with where in the sequence (in that PDB code).
(I wrote this in the reply box, so there may be a typo or two).
PyMOL has a selection algebra called pepseq, which matches missing atoms, but I was a whimp and played it safe by being more verbose. Also, a for loop can be used instead of an iterate. But I think it is slower.

atoms = pymol.cmd.get_model(f"pepseq {target} and name CA")
for atom in atoms.atom:
      if atom.ss == ' ...

Uniprot

In Uniprot, if a protein is crystallised, there will be the feature entries corresponding to the secondary structure. It is sheet or helix. No turns, pi-helices or 3_10 helices. There is no Biopython method for parsing Uniprot XML. I wrote myself a parser and it is not too horrid thing to do —although I monkeypatched elementtree. Alternatively, xmlschema.XMLSchema.to_dict can be used to convert the XML to a dictionary. The advantage of Uniprot is that the residues ids will not require correcting. The disadvantage is parsing XML.

My parser is specific for my application and it is poorly documented (so would not recommend it) but but... I can share it and human data. I'd be used like:

 import os
 from michelanglo_protein import ProteinCore
 ProteinCore.settings.startup(data_folder='protein-data')
 p = ProteinCore().gload(file=os.path.join('gpickle','Q86V25.gpz'))
 #print(p.asdict())
 if 'helix' in p.features:
        n = p.sequence.find( ...
        for ... p.features['helix'] #('helix', 'turn', 'strand')
        if x < n < y ...


Cannot blast database created with blastdb_aliastool: BLAST Database error: No alias or index file found

The problem was that the database must be specified without suffix, i.e. blastn -db testdb as opposed to blastn -db testdb.nal.

By the way, it produces exactly the same output as

blastn -db "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -query query_single_nucl.fasta -outfmt 6

and the performance is not better.


isSpike function in SingleCellExperiment package is deprecated?

isSpike was provided by SingleCellExperiment. It is now deprecated. The notes in SingleCellExperiment indicate that one should use altExps

I am struggling with a similar problem. Can't compile a package (powsimR) because a dependency uses isSpike (not used in the code directly), so I have to hunt down where it is used by one of a very long list of dependencies.


Huge variability in blastn speed between same short sequence job

First of all, if you're using an online service like NCBI's BLAST, then the time to completion will depend on the load the service is experiencing. Your job is ut in a queue and will be started when it's its turn. If the service is under heavy load, then its turn might come later than you expect. Then, if the service is under heavy load, even after launching your job, it might take longer to finish because it won't have enough resources to run quickly.

Since you say that the same job finished quickly on a new tab, it sounds like either the first job was last in the queue or there was just a bug that caused it not be launched or something like that. We can't really help with things that resolve themselves like that.


Database of copy number alterations in different cancer types

I don't know if there's a database that does exactly what you want, but there are some places that might help you figure this out, especially if you already have a list of CNAs/genes/regions in mind.

ICGC is the first that comes to mind, as it has samples from many types of cancer and copy number data for most of them. Depending on what you want though, you might have to apply for access to download the data you want. It's rather annoying, and their data download tools are atrocious to use, but they have a lot of data.

DGV might be another option, though their data browser is rather dated feeling. They publish a lot of gold-standard germline CNVs and such.

COSMIC is a large database for cancer-associated mutations, including CNVs, so it might be worth a look for you.


Publicly available, free, complete database for antibiotics names and classes?

ChEBI

It should be possible to get this information from the ChEBI database, see the exported tables. You could download the ontologies (in OWL / OBO), parse them using some ad-hoc parser or using a dedicated library (e.g., Pronto), build a directed acyclic graph based on the is_a relationship, and extract the antibiotic subtree.

KEGG

You could also use KEGG: Kyoto Encyclopedia of Genes and Genomes, but a subscription is probably required. The exported drug table should be easy to parse (e.g., in Python).

For a quick verification of the completeness, you can download the table, filter out everything except antibacterial agents and check the number of records:

curl "ftp://ftp.genome.jp/pub/kegg/medicus/drug/drug" > keg_all.txt

cat keg_all.txt \
   | perl -pe 's/\n/BREAK/g' \
   | perl -pe 's@///BREAK@///\n@g' \
   | grep 'CLASS       Antibacterial' \
   | perl -pe 's/BREAK/\n/g' \
   > keg_antibacterial.txt

cat keg_antibacterial.txt | grep "///" | wc -l

There are 646 records. However, some of them are almost identical (e.g., only 1 atom difference).

If you add | grep -i "antibiotic" \ to the pipeline, then you get only 445 records.


Mapping Reads to Known Gene Paralogs with Long Read Technology

Question 1

This command will get you all of supplementary alignments for the reads. This isn't exactly what you want though. You want all of the reads that have more than one mapping.

samtools view -f 2048 -h myfile_sorted.bam > supp_only.bam

This bash script returns a bam file that contains all alignments for reads that have multiple mappings. It works by only selecting reads that have supplementary alignments (2048 flag), or are not the primary alignment (256 flag), then using picard tools to filter the bam file.  Using picard tools is the fastest way I have seen to filter a bam file given a list of reads.

#/bin/bash
BAMFILE=/path/to/my_file.bam
OUTFILE=/path/to/multiply_mapped.bam
PICARD_JAR=/path/to/picard.jar
TEMP=multiply_mapped.txt

for filt_int in 256 2048; \
do samtools view -f "$filtint" ${BAMFILE} | \
  grep -v '#' | cut -f1; done | sort | uniq > ${TEMP}; \
java -jar ${PICARD_JAR} FilterSamReads \
  I=${BAMFILE} \
  O=${OUTFILE} \
  READ_LIST_FILE=${TEMP} \
  FILTER=includeReadList
rm ${TEMP}

This command took 21 seconds to execute for my bam file with ~300,000 long reads mapped and several repetitive regions.

Question 2

Filtering the bam file should be pretty standard, as addressed above in my response to Q1. However, different long-read aligners are more sensitive to distantly-related paralogs. You might consider trying the blasr aligner since it is tuned for PacBio long reads, or try bwa-mem to compare your results with minimap2. Realistically though, most of these aligners will probably behave similarly to one another.

The fact that there are paralogs does not matter to the aligners. The aligner just sees these regions as two different places to which a single read can map just about equally well. When there are multiple similar regions to which a read can map, as in the case of paralogs, the aligner will report both alignments. I do not know of any long-read paralog-aware aligners.

As to why none of the reads mapped to the paralogs: Perhaps they have low mappability, as in they are repetitive? Also, maybe the capture didn't work properly? It looks like there is serious sequencing bias to the dataset given the uneven coverage.

Question 3

Are the reads that are mapping actually 4kb, or are the alignments only 4kb long? You may have many 6kb-long sequences that have been clipped at either the 5' or 3' end. If this is not a heavily-curated genome then the assembly may be incorrect around the paralogs, preventing proper mapping.

One suggestion is to take whole-genome long reads from the same organism and to map them. Then check if the mapping to this region is equally strange compared to your capture data. This will give you a better idea of what caused the uneven coverage and seeming lack of coverage around the paralogs.


Where can I get the population allele frequency vcf file?

On the GATK forum they've recommended the population stratified VCF file for this purpose.


In R:

flist = list.files("*.txt", path=".")
fl = lapply(flist, function(x) read.delim(x, header=F))
m = sapply(fl, function(x) x$V2)
row.names(m) = fl[[1]]$V1
colnames(m) = flist

That will give you an appropriately formatted matrix named m in R, which is probably where you want to continue processing things anyway.


how to merge two fasta based on full fasta header and remove the duplicate

You can use the FastaToTbl and TblToFasta scripts I have posted here for this:

FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa

If you run this on your example, you will get:

$ FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa
$ cat merged.fa
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATG
GTGCGGGAACGCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGT
TTGGATGTCTCTGGCGACGG
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGC
GAACGCTACGGGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGAT
GTCACCGACGACGGTGGGTG
>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAG
GTCGCGGTCGGCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTC
GACCGCCTCGAAGCCCGCCT
>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAG
GTCGCGGTCGGCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTC
GACCGCCTCGAAGCCAGCGC
>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGC
CAGTTGGTACGGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGAC
GTCGGCCTCGATGTCACGGG

This will remove any cases where both the header and the sequence are identical. If you instead want to do this based on the header only, not taking the sequences into account, you can just process both files and only print sequences from the 2nd if their ID line hasn't been seen before. For example, with awk:

$ awk '{
        if(NR==FNR){
            if(/^>/){seen[$0]++;}
            print
        }
        else{
            if(/^>/){ skip=seen[$0] }
            if(!skip){ print }
        }
       }' file1.fa file2.fa
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT


I found the way to download the genes and pathways thanks to this answer using KEGGREST.

It seems that the API allows to fully download the database because it allows to download info, list, find, get, conv, link, ddi.

In fact in the message of 2011-05-21 it mentions that the API includes more options to download the KGML files that previously where only possible to get via the FTP site.


Do variant calls change when you call from CRAM?

By default, a CRAM you create with samtools is lossless. It typically halves the input BAM in terms of file size. If you want to compress more, you can let samtools convert most read names to integers. You won't be able to tell optical duplicates from read names, but this is a minor concern. You can also drop useless tags depending on your mapper and the downstream caller in use. For cancer data, I wouldn't reduce the resolution of base quality without comprehensive benchmarks. Unfortunately, base quality takes most of space in CRAM. Discarding the original read names and some tags probably won't save you much space.


Calculating bit score: How do you find lambda and K?

NCBI BLAST uses tables of pre-computed values for lambda/K for gapped protein and nucleotide comparison and their values depend on the amino acid substitution matrix or the substitution scores + gap penalties.

You can find these tables (from line 589 for blastn) in the file c++/src/algo/blast/core/blast_stat.c from the blast source package.


a package with randomization/permutation tests for genomic analyses

I would simply look at co-variance across the entire dataset and examine values >0.8 in the first instance rather than select 1000 genes.

With your approach sampling 1000 genes can be criticised for selective sampling and the way to get around this is to perform replicates, with (bootstraps) or without replacement. You could produce e.g. 1000 replicates of 1000 genes and perform co-variance on each replicate to examine Y = X.

Python has a specific permutation/combinations module call itertools, which is different to random resampling and doesn't need replicates

from itertools import combinations
X = [1,2,3]
Y = [1,2,3]
print (list(combinations(iterable, 1000) + X + Y))

The standard way is to sample is,

from random import sample
mygenes = [1, 2, 3, 4, 5]

print(sample(mygenes,3))

Output

[2, 4, 5]

There is nuance and further funcationality ... I cannot remember is "sample" is with or without replacement (its quite a big deal)


Visualization tools for 3C/Hi-C long-range interaction data?

I would also recommend two very recent Hi-C visualization frameworks (with some public data available in both): HiGlass and JuiceBox.


A Database for Plant Metabolome or Secondary Metabolites?

PhytoMine is an InterMine interface to data from Phytozome, the Plant Comparative Genomics portal of the Department of Energy's Joint Genome Institute, which hosts 81 assembled and annotation genomes, from 78 viridiplantae species.


Finding common and unique data set by comparing two files based on their column and to split the columns multiple strings to print in output

One way to do this is to use a dictionary. You can read the first file line by line, and for each line, use "chr:position" as a key, and the corresponding line as the value. If you want to handle duplicates, you can make the value a list and store multiple lines with duplicated "chr:position" as elements in that list.

After constructing this dictionary, read the second file line by line. For each line, set "chr:position" as the so-called name of that line. Then using this "chr:position", you can easily find the corresponding line in the dictionary. Append the line from the second file to that line from the first file.

After finishing reading the second file, you are ready to write the new file to output. To do so, iterate over the keys in the dictionary, and write the value to each key to a new file. Remember that each value in that dictionary is now a new line.

I choose to use a dictionary because they are basically a hash table, meaning that using the key you can instantly read that value in the dictionary, instead of iterating every element in a list to find that element. This is very handy when you have a large dictionary.


It would be useful if you show us some line of your CADD.bed and Input.vcf and how your desired output should look like.

bedtools is complaining about non positional records in CADD.bed, which means there are lines in this file where the first three columns aren't CHROM, START, END.

I guess you will have some header lines in the beginning that aren't introduced by an #. Just delete them (or add in # to the beginning of the line) and try again.

fin swimmer


bigWig to bed for regions above/below threshold

You can do that with a bit of python:

#!/usr/bin/env python
import pyBigWig
threshold = 10  # Change me
bw = pyBigWig.open("file.bw")  # Change me
of = open("regions.bed", "w")  # Change me

for chrom, len in bw.chroms().items():
    intervals = bw.intervals(chrom)
    for interval in interals:
        if abs(interval[2]) > threshold:
            of.write("{}\t{}\t{}\n".format(chrom, interval[0], interval[1]))
bw.close()
of.close()

You'll need to install pyBigWig (it's available via pip or conda install) and change the lines with Change me in them.


Which tools can detect chimeric RNA (fusion genes) from WGS or RNA-Seq data?

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

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

Other useful programs:

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

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


How to run HYPHY on multiple files non-interactively?

TL;DR

I made the following bash script to generate an input file for HYPHY.

#!/bin/bash
cat << EOF
inputRedirect = {};
inputRedirect["01"]="Universal"; // genetic code
inputRedirect["02"]="$(readlink -f $1)"; // codon data
inputRedirect["03"]="$(readlink -f $2)"; // tree
inputRedirect["04"]="${3:-All}"; // Test for selection on a branch
inputRedirect["05"]=""; // complete selection

ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles/BUSTED.bf", inputRedirect);
EOF

If you want to test all the branches jointly run:

./gen_busted.sh alignment.phy tree.nwk > script.bf

Or for a particular branch:

./gen_busted.sh alignment.phy tree.nwk branchLabel > script.bf

And then run:

HYPHYMP script.bf

Now you can use the script to generate a .bf file for every sequence alignment you have, and run HYPHY for every file.

Longer version

I will describe how to make a similar script for any model, such as MEME, RELAX or PARRIS.

First, find the .bf script which performs the analysis. On GNU/Linux it should be at HYPHY_INSTALLATION_PATH/lib/hyphy/TemplateBatchFiles (it might be different between different OS). In case of BUSTED your script has a name BUSTED.bf.

Now run this model in the interactive mode in order to record all the inputs required by the model. Make sure to specify full paths to all the files.

HYPHYMP HYPHY_INSTALLATION_PATH/lib/hyphy/TemplateBatchFiles/BUSTED.bf

In case of BUSTED the inputs are:

Genetic code.
Codon data (sequence alignment in phylip format).
Tree (in newick format).
Which branches to test.

Now for every input alignment you need to generate a batch file which looks like this:

inputRedirect = {};
inputRedirect["01"]="Universal"; // genetic code
inputRedirect["02"]="/path/to/alignment.phy"; // codon data
inputRedirect["03"]="/path/to/tree.nwk"; // tree
inputRedirect["04"]="All"; // Test for selection on all branches
inputRedirect["05"]="BRANCH1"; // Test for selection on branch1
inputRedirect["06"]="BRANCH2"; // Test for selection on branch2
inputRedirect["07"]=""; // complete selection

ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles/BUSTED.bf", inputRedirect);

Now you can create a script which generates the file. The most important thing, don't forget to use full paths in the .bf file.


Program to make a haplotype network for a specific gene

This is a haplotype map, where each node is proportional to the frequency of a given allele. The formal tool to do this is Network 5 found here. The problem with manaul haplotype maps is the area of the node should be proportional to the frequency, if this is not accurate it could be perceived to be misleading.

The approach for looking at each parsimonious state along a branch is a phylogenetics method which is performed via MacClade 4.0 ... requires an old Mac though.


What is the best way to get a large number of RNA seq data from SRA in Python without being denied access

I once needed to download more than 1000 files from SRA, and I found a blog post that is very useful (https://reneshbedre.github.io/blog/fqutil.html). The idea is to use aspera to fetch SRA files. This is fast becasue of aspera. Then use fasterq-dump to convert SRA file to fastq. fasterq-dump is faster than fastq-dump partly because it does not compress the output file. I think compressing is the rate-limiting step here. After this, you can use pigz, which is a multi-threading version of gzip to compress the fastq files. Note that if you want to download files to a computing cluster, you can submit jobs for fasterq-dump and pigz. This makes things much faster. You cannot submit jobs for fetching SRA, because you need to connect to internet.


How to know if ID correspond to nucleotide sequence or protein with a script

If these are all GenBank or RefSeq accessions, you can use Entrez Direct for this as shown below:
$ cat accs.txt
ATO98108.1
ATO98120.1
ATO98132.1
AVP78031.1
AVP78042.1
$ cat accs.txt | epost -db nuccore | efetch -format acc
## no output because none of them are nucleotide accessions
$ cat accs.txt | epost -db protein -format acc | efetch -format acc
AVP78042.1
AVP78031.1
ATO98132.1
ATO98120.1
ATO98108.1

NOTE: This will work only if the accessions are currently live because epost does not find any suppressed accessions. For example:
$ cat accs.txt
NM_002826.3
NM_002826.4
NM_002826.5
$ cat accs.txt | epost -db nuccore -format acc | efetch -format acc
NM_002826.5

Here, all three accessions are valid nucleotide accessions but only the last one, NM_002826.5, is alive.
An alternate way is to use the accession prefixes defined here and come up with an appropriate regular expression query.


How to select high quality structures from the Protein Data Bank?

There is a very nice database, pdbcull (also known as the PISCES server in the literature). It filters the PDB for high resolution and reduced sequence identity. It also seems to be updated regularly. Depending on the cut-offs, you get between 3000 and 35000 structures.

If you are specifically interested in rotamers, you may want to look at top8000 instead, where they have checked for high resolution, and good MolProbity scores. They also provide a rotamer database.

PDB also provides their own clustering. They first cluster the sequences, and then extract a representative structure for each one, based on the quality factor (1/resolution - R_value). This has the advantage of being comprehensive, but you will have bad structures when no good ones were ever obtained.


scanBam from Rsamtools is not importing one of my reads into R

Note the flag; that read is mapped in a reverse-complemented manner, so isMinusStrand = FALSE is filtering it out. I tested this by making a BAM file with only that read:

@SQ SN:I    LN:1000000
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   00  *   *

Then in R:

> library(Rsamtools)
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = TRUE)))[[1]]
> blah$qname
[1] "000000F|arrow"
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = FALSE)))[[1]]
> blah$qname
character(0)


How to find correlation between two specific genes in same dataset?

I suppose in an idealized world you'd use something like salmon to get TPMs, which can account for differences in effective gene length by sample. However that's quite likely to be overkill.
Not necessary, you're using a with-sample ratio, which wouldn't be changed by normalization.

Here is a small example in R, assuming your matrix or data frame is named m:

first = which(row.names(m) == "ESCO1")
second = which(row.names(m) == "ESCO2")
plot(m[first,], m[second,])

or cor(m[first, ], m[second,]). You'll want a rank-based correlation.


Using external list of PCs for clustering

There is a way to do this, and even better--there is documentation for how to do it! No surprise coming from the Satija Lab. In the vignette they perform multidimensional scaling, but the idea is the same. cmdscale() returns the cell embeddings. SetDimReduction() is the Seurat function you are looking for. No manual editing using @ required.

The authors use add mds embeddings, but if you want to add your own pca embeddings swap out "mds" for "pca" and mds for URD_object@pca.load.

# Before running MDS, we first calculate a distance matrix between all pairs
# of cells.  Here we use a simple euclidean distance metric on all genes,
# using object@scale.data as input
d <- dist(x = t(x = pbmc@scale.data))
# Run the MDS procedure, k determines the number of dimensions
mds <- cmdscale(d = d, k = 2)
# cmdscale returns the cell embeddings, we first label the columns to ensure
# downstream consistency
colnames(x = mds) <- paste0("MDS", 1:2)
# We will now store this as a new dimensional reduction called 'mds'
pbmc <- SetDimReduction(object = pbmc, reduction.type = "mds", slot = "cell.embeddings",
    new.data = mds)
pbmc <- SetDimReduction(object = pbmc, reduction.type = "mds", slot = "key",
    new.data = "MDS")

See the dimensional reduction vignette for more.


I know that structure way too well —SARS-COV-2 MPro.
In PyMOL command line
 remove resn 02J or 010

or
 cmd.remove('resn 02J or 010')

In python
 import pymol2
 with pymol2.PyMOL() as pymol:
      pymol.cmd.fetch('6LU7')
      pymol.cmd.remove('resn 02J or 010')
      pymol.cmd.save('6LU7.mod.pdb')
      #or pdbblock = pymol.cmd.get_pdbstr()

From the sounds of it you'd like to make the peptidomimetic a peptide. If you wanted to do the horrid thing of making PJE, the capped side chain glutamine ketamide, into glutamine, the command to rename atoms is a multitude of cmd.alter('resn PJE and name <atomname>', 'name="<newname>") followed by a cmd.alter('resn PJE', 'resn="GLN"') and lastly cmd.sort() to commit the changes. Not pretty, but is all in PyMOL.


Parsing SAM/BAM files for additional information

Other options... samtools stats maybe gives you some of the stats you are after. Chances are you already have samtools so nothing additional to install. Also useful may be pysamstats, in particular the variation sub-command.


It doesn't look like this data has been clustered, and you are plotting the 'value', rather than the cluster. If you want to represent values on a tSNE using Seurat, a better function to use is the FeaturePlot, which will colour the points in a more meaningful colour gradient based on their value.
[It's possible this function wasn't available at the time this question was written]
If, instead, you want to plot clusters, then you should carry out some clustering on the data first using FindNeighbors/FindClusters.
There's a great introductory tutorial to guide you into using Seurat on the Satija lab page:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html


Calculate genome coverage and depth from alignment

You could use samtools coverage as explained in the manual of samtoools.
Here is a example which is also described on the manual site.
samtools coverage -r chr1:1M-12M input.bam

#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8


differential expression multiple studies

Adding experiments 1 & 2 should be fine without modification on the computer side of things, as long as your sample preparation is consistent. To be ultra-cautious, you could add in the experiment as an additional covariate for the differential expression model.

I'm less confident about adding in TCGA data; it's probably going to look better than your own samples, which might cause problems in calculating differential expression from raw reads. You can do a separate normalisation/correction of the TCGA data, but that will limit the inferences that can be made about differential expression -- without raw count information, the shot noise associated with read sampling cannot be included in the statistical model. I notice that TCGA has already done their own filtering / expression, and that would be preferable to use (mRNA sequencing, data subtype 'Gene'):

https://tcga-data.nci.nih.gov/docs/publications/tcga/datatype.html

However, even with data that has been properly normalised and filtered to exclude batch effects, there will still be population differences that influence results. My recommendation would be to do some subsampling to reduce the population effects, something like at least the following comparisons:

6 Vs 100 + 5 normals
6 Vs 100 normals
6 Vs 50 + 3 normals [randomly sampled from the two normal groups]

And only consider differential expression that is consistently observed in all three of these comparisons. Note that the third "subsampled" dataset can be generated more than once, taking different subgroups each time.

If you are going to use separately normalised data, then DESeq2 may not be appropiate (because it expects count-level data). Something like Limma could work for DE analysis, see chapter 15 of the user's guide for RNASeq analysis. Unfortunately I'm not very familiar with using Limma, so can't provide any advice on things like what to watch out for in the input data and/or results.


How to remove batch effect from TCGA and GTEx data

This has been asked and answered here and Biostars. My take home message from many similar posts is "not to bother with correcting batch effects when performing differential expression analysis but to include 'batch' as a covariate". The latter accounts for the batch effects.
EDIT: This does not really answer OP's question, as @swbarnes2 and @ATpoint pointed out, they have a completely confounded design and the linear regression model would not accept such a model matrix.


Error while calling bcftools mpileup - Failed to open -: unknown file type

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

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

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


Human Herpes Virus, fusion motif in protein vs. SARS-CoV-2?

I can answer that ... I'm afraid no. What you are asking is whether there is any similarity in their cellular receptor binding or "fusion loops" because the motifs look a bit similar.
These are very different viruses and you need to appreciate just how different they are,
Human herpes is a double stranded DNA virus resulting in a chronic infection. It has a vast genome of >150kb, which includes cytokine mimics and proteins that we have no idea what they do. The protein you are referring to is likely a surface antigen in a virus covered in different surface antigens. However, this protein is involved in blocking the beta inferon response is considered a protease - not fusion loop and appears to have some sort of complex replication function involving the capsid. Finally, there is also the question of which human herpes 1 to 8. Here you are referring to HHV1 but for example HHV3 is chickenpox (okay its an acute not chronic infection) and HHV5 is cytomegliovirus (congenitcal infection).
COVID-19 is a simple virus of a 30Kb genome of a single stranded positive sense RNA, causing an acute infection. I agree the Spike protein is much bigger than most RNA viruses, but the receptor-binding is tightly defined with ACE-2. It forms a classic RNA icosahedral virus, singly covered in Spike protein. I have no doubt it will go through conformational change like other RNA viruses in receptor binding but the protein involved is singularly involved in cell entry - thats its job.
Summary
There are large differences between HHV1 and COVID-19 for,

Viruses; DNA double stranded vs RNA single stranded
Infections; chronic vs. acute
Transmission; vertical/integration? direct contact vs. respiratory
Ultimately very different proteins
Finally very different functions

At best this would be heavily  scrutinised as a virological assumption.
Bioformatics If you really wanted to follow this up you could compare the structures around these motifs, because COVID-19 is close to SARS and assess this possibility ... but you would need to be careful about the conformational changes Spike goes through in receptor-binding.
Advice before investing large amounts of time in a bioinformatics question, it is well-worth understanding the biological rationale for the question. In this instance that rationale appears to be lacking. If you found a comparison for example with gamma-coronaviruses, which was absent from many beta-coronaviruses that would be interesting. Here there is simply too much possibility of a chance similarity and given the size HHV1 that cannot be easily overlooked. However it is good that you are thinking and questioning, if you apply the same techniques, e.g. psi-blast, protein-blast into different questions ... that can work.


Genome assembly of SRR12196449 with SPAdes

Update 2:
I looked into this a little more, with the various data sources.
This is related in part to the answer submitted by OP juanjo75es, in addition to discussion on chat. I don't entirely understand the logic, but the general thrust seems to be that SPAdes makes weird assemblies somehow.
Some notes that I made:
REFERENCE ASSEMBLIES

FIV sequence U11820.1 was deposited in 1996, before SPAdes existed. Unclear assembly method.
FIV sequence MN630242 was deposited in 2020, but was assembled with CLC workbench.
These two sequences align together ok with minimap2. Low exact identity, it's true. This is the PAF output:

MN630242.1      8977    194     8964    +       U11820.1        9462    575     9326    1130    8778    60      tp:A:P  cm:i:116        s1:i:1125
       s2:i:55 dv:f:0.1761

These two assemblies are syntenic, and both have similar gene calls using prokka. GFFs:

MN630242.1      Prodigal:2.6    CDS     253     1605    .       +       0       ID=AMGANPBD_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00001;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     1656    4868    .       +       0       ID=AMGANPBD_00002;eC_number=3.6.1.23;Name=dut;gene=dut;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:Q2YRG4;locus_tag=AMGANPBD_00002;product=Deoxyuridine 5'-triphosphate nucleotidohydrolase
MN630242.1      Prodigal:2.6    CDS     4861    5616    .       +       0       ID=AMGANPBD_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00003;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     5891    8461    .       +       0       ID=AMGANPBD_00004;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00004;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     8626    8790    .       +       0       ID=AMGANPBD_00005;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00005;product=hypothetical protein

U11820.1        Prodigal:2.6    CDS     634     1983    .       +       0       ID=EANIPDKN_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00001;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     1995    5246    .       +       0       ID=EANIPDKN_00002;eC_number=3.6.1.23;Name=dut;gene=dut;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:Q2YRG4;locus_tag=EANIPDKN_00002;product=Deoxyuridine 5'-triphosphate nucleotidohydrolase
U11820.1        Prodigal:2.6    CDS     5239    5994    .       +       0       ID=EANIPDKN_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00003;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     6269    8830    .       +       0       ID=EANIPDKN_00004;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00004;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     8904    9152    .       +       0       ID=EANIPDKN_00005;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00005;product=hypothetical protein

I can share FAA files of the proteins if needed.
READ ASSEMBLY
I also assembled the indicated reads using SPAdes. For reference it is a ~9Kbp virus, but this is a 4.3Mbp assembly. There is a lot of non-virus sequence in there. The second-largest contig is a shuttle vector. the third-largest contig is cat (host). Many more are cat, so I think it's pretty cat-oriented. The original authors used CLC workbench, so I guess that just worked a lot better in this instance, even in the presence of all the contamination. Unclear why, it appears that CLC works similarly to OP's assembly tool.
I aligned these to MN630242.1 with minimap2. If I understand OP, they are unhappy about the mapping of these contigs to this genome reference. The identities are mostly high (though there is indeed some oddness with the largest contig, which only finds very low coverage- at the same time that there are overlapping contigs with very high ID?), and they cover the entire reference genome:
MN630242.1  8977    2155    3199    -   NODE_14_length_1054_cov_3786.620280 1054    2   1046    1021    1044    60  tp:A:P  cm:i:189    s1:i:1021   s2:i:378    dv:f:0.0031
MN630242.1  8977    7893    8626    -   NODE_25_length_743_cov_3814.258117  743 8   739 713 733 60  tp:A:P  cm:i:129    s1:i:713    s2:i:318    dv:f:0.0045
MN630242.1  8977    72  4904    -   NODE_1_length_4942_cov_25.814123    4942    34  4863    628 4832    41  tp:A:P  cm:i:7s1:i:628  s2:i:512    dv:f:0.1703
MN630242.1  8977    3376    3897    +   NODE_39_length_526_cov_5254.155388  526 3   524 512 521 0   tp:A:S  cm:i:8s1:i:512  dv:f:0.0023
MN630242.1  8977    6270    6673    +   NODE_82_length_409_cov_1.744681 409 3   406 357 403 0   tp:A:P  cm:i:53 s1:i:357    s2:i:351    dv:f:0.0185
MN630242.1  8977    6852    7223    -   NODE_87_length_381_cov_23.334646    381 6   376 355 371 2   tp:A:P  cm:i:5s1:i:355  s2:i:351    dv:f:0.0079
MN630242.1  8977    6187    6558    -   NODE_88_length_380_cov_2338.675889  380 7   378 351 371 0   tp:A:S  cm:i:5s1:i:351  dv:f:0.0085
MN630242.1  8977    6852    7308    +   NODE_56_length_471_cov_0.933140 471 15  471 351 456 0   tp:A:S  cm:i:43 s1:i:351    dv:f:0.0397
MN630242.1  8977    7332    7707    +   NODE_86_length_389_cov_3506.202290  389 5   380 334 375 17  tp:A:P  cm:i:5s1:i:334  s2:i:309    dv:f:0.0111
MN630242.1  8977    5181    5568    -   NODE_60_length_464_cov_1.005935 464 50  437 314 387 17  tp:A:P  cm:i:41 s1:i:314    s2:i:282    dv:f:0.0357
MN630242.1  8977    7277    7707    +   NODE_76_length_438_cov_0.919614 438 5   435 309 430 0   tp:A:S  cm:i:40 s1:i:309    dv:f:0.0437
MN630242.1  8977    6953    7285    +   NODE_94_length_345_cov_2.903670 345 8   340 301 332 0   tp:A:S  cm:i:50 s1:i:301    dv:f:0.0110
MN630242.1  8977    5744    6056    -   NODE_113_length_316_cov_2.169312    316 4   316 289 312 35  tp:A:P  cm:i:4s1:i:289  s2:i:244    dv:f:0.0096
MN630242.1  8977    5615    5927    +   NODE_108_length_322_cov_1103.400000 322 6   318 283 312 58  tp:A:P  cm:i:4s1:i:283  s2:i:209    dv:f:0.0143
MN630242.1  8977    6543    6813    -   NODE_147_length_279_cov_2592.519737 279 8   278 261 270 21  tp:A:P  cm:i:4s1:i:261  s2:i:236    dv:f:0.0032
MN630242.1  8977    8544    8872    +   NODE_100_length_333_cov_2259.189320 333 1   329 252 328 0   tp:A:P  cm:i:3s1:i:252  s2:i:245    dv:f:0.0318
MN630242.1  8977    7601    7853    -   NODE_200_length_255_cov_25.125000   255 0   252 250 252 12  tp:A:P  cm:i:3s1:i:250  s2:i:236    dv:f:0.0017
MN630242.1  8977    6648    6968    -   NODE_106_length_324_cov_2.538071    324 4   324 248 320 0   tp:A:P  cm:i:3s1:i:248  s2:i:248    dv:f:0.0303
MN630242.1  8977    5940    6193    -   NODE_162_length_266_cov_1809.287770 266 6   259 246 253 5   tp:A:P  cm:i:4s1:i:246  s2:i:240    dv:f:0.0033
MN630242.1  8977    5416    5664    -   NODE_264_length_255_cov_9.531250    255 5   253 243 248 0   tp:A:P  cm:i:4s1:i:243  s2:i:242    dv:f:0.0033
MN630242.1  8977    7221    7466    -   NODE_247_length_255_cov_12.539062   255 1   246 241 245 0   tp:A:P  cm:i:4s1:i:241  s2:i:241    dv:f:0.0032
MN630242.1  8977    4830    5071    -   NODE_1121_length_245_cov_11.711864  245 4   245 239 241 0   tp:A:P  cm:i:4s1:i:239  s2:i:235    dv:f:0.0042
MN630242.1  8977    7769    8013    +   NODE_677_length_251_cov_3.951613    251 6   249 238 244 7   tp:A:P  cm:i:4s1:i:238  s2:i:230    dv:f:0.0075
MN630242.1  8977    4901    5148    +   NODE_270_length_255_cov_8.726562    255 8   255 238 247 0   tp:A:P  cm:i:4s1:i:238  s2:i:234    dv:f:0.0073
MN630242.1  8977    8705    8948    -   NODE_382_length_254_cov_5.763780    254 7   250 234 243 15  tp:A:P  cm:i:4s1:i:234  s2:i:217    dv:f:0.0062
MN630242.1  8977    6065    6307    +   NODE_410_length_254_cov_3.488189    254 1   243 227 242 10  tp:A:P  cm:i:3s1:i:227  s2:i:210    dv:f:0.0122
MN630242.1  8977    5062    5239    -   NODE_7003_length_180_cov_9.207547   180 3   180 172 177 17  tp:A:P  cm:i:2s1:i:172  s2:i:157    dv:f:0.0066
MN630242.1  8977    11  136 +   NODE_21321_length_128_cov_1195.000000   128 2   127 125 125 3   tp:A:P  cm:i:2s1:i:125  s2:i:123    dv:f:0

Only a (relatively) few contigs align well (figure), as might be expected with heavy contamination, and of those some are overlapping:

To look into the odd behavior of the biggest contig, I ran prokka on it too and found the expected genes for its position, especially the largest gene in the virus, dut (Deoxyuridine 5'-triphosphate nucleotidohydrolase). I then took the protein sequences of the 3 assemblies and aligned them with clustalo. They're all quite similar:
CLUSTAL O(1.2.4) multiple sequence alignment

NODE1_KJHFFCBH_00001         ------MEKRPEIQIFVNGHPIKFLLDTGADITILNRKEFIIGNSIENGKQNMIGVGGGK
MN630242_AMGANPBD_00002      -------------MIFVNGYPIKFLLDTGADITILNRRDFQVKNSIENGRQNMIGVGGGK
U11820.1_EANIPDKN_00002      MGTTTTLERRLEIQIFVNGHPIKFLLDTGADITILNRKDFQIGNSIENGKQNMIGVGGGK
                                           *****:*****************::* : ******:**********

NODE1_KJHFFCBH_00001         RGTNYINVHLEIRDENYKTQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISEKI
MN630242_AMGANPBD_00002      RGTNYTNVHLEIRDENYKTQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISDKI
U11820.1_EANIPDKN_00002      RGTNYINVHLEIRDENYRMQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISEKI
                             ***** ***********: **************************************:**

NODE1_KJHFFCBH_00001         PIVKVRMKDPTQGPQVKQWPLSNEKIEALTEIVERLEQEGKVKRADPNNPWNTPVFAIKK
MN630242_AMGANPBD_00002      PIVKVKMKDPNKGPQIKQWPLSNEKIEALTEIVERLEKEGKVKRADPNNPWNTPVFAIKK
U11820.1_EANIPDKN_00002      PIVKVRMRDPIQGPQVKQWPLSNEKIEALTDIVERLESEGKVKRADPNNPWNTPVFAIKK
                             *****:*:** :***:**************:******.**********************

NODE1_KJHFFCBH_00001         KSGKWRMLIDFRVLNKLTDKGAEVQLGLPHPAGLQWKKQVTVLDIGDAYFTIPLDPDYAP
MN630242_AMGANPBD_00002      KSGKWRMLIDFRELNKLTEKGAEVQLGLPHPAGLQMKKQVTVLDIGDAYFTIPLDPDYAP
U11820.1_EANIPDKN_00002      KSGKWRMLIDFRVLNKLTDKGAEVQLGLPHPAGLQMKKQVTVLDIGDAYFTIPLDPDYAP
                             ************ *****:**************** ************************

NODE1_KJHFFCBH_00001         YTAFTLPRKNNAGPGKRYVWCSLPQGWVLSPLIYQSTLDNILQPYIKQNPELDIYQYMDD
MN630242_AMGANPBD_00002      YTAFTLPRKNNAGPGRRYIWCSLPQGWILSPLIYQSTLDNIIQPFIRQNPQLDIYQYMDD
U11820.1_EANIPDKN_00002      YTAFTLPRKNNAGPGRRYVWCSLPQGWVLSPLIYQSTLDNILQPFIKQNSELDIYQYMDD
                             ***************:**:********:*************:**:*:** :*********

NODE1_KJHFFCBH_00001         IYIGSNLSRQEHKQKVEELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWSIQQKQL
MN630242_AMGANPBD_00002      IYIGSNLNKKEHKEKVGELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWTIQQKQL
U11820.1_EANIPDKN_00002      IYIGSNLSKKEHKQKVEELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWSIQQKQL
                             *******.::***:** ************************************:******

NODE1_KJHFFCBH_00001         EIPERPTLNDLQKLAGKINWASQTIPKLSIKALTHMMRGDQKLDSIREWTEEAKKEVQKA
MN630242_AMGANPBD_00002      DIPEQPTLNELQKLAGKINWASQAIPDLSIKALTNMMRGNQNLNSIREWTKEARLEVQKA
U11820.1_EANIPDKN_00002      EIPERPTLNELQKLAGKINWASQTIPDLSIKELTNMMRGDQKLDSIREWTVEAKREVQKA
                             :***:****:*************:**.**** **:****:*:*:****** **: *****

NODE1_KJHFFCBH_00001         KEAIEKQAQLYYYDPNRELYAKISLVGPHQLCYQVYHKNPEQILWYGKMNRQKKRAENTC
MN630242_AMGANPBD_00002      KKAIEGQVQLGYYDPSKELYAKLSLVGPHQISYQVYQRNPERILWYGKMSRQKKKAENTC
U11820.1_EANIPDKN_00002      KEAIEKQAQLNYYDPNRGLYAKLSLVGPHQICYQVYQKNPEHILWYGKMNRQKKKAENTC
                             *:*** *.** ****.: ****:*******:.****::***:*******.****:*****

NODE1_KJHFFCBH_00001         DIALRACYKIREESIVRIGKEPVYEIPASREAWESNIIRSPYLKASPPEVEFIHAALNIK
MN630242_AMGANPBD_00002      DIALRACYKIREESIIRIGKEPKYEIPTSREAWESNLINSPYLKAPPPEVEYIHAALNIK
U11820.1_EANIPDKN_00002      DIALRACYKIREESIIRIGKEPMYEIPASREAWESNLIRSPYLKAPPPEVEFIHAALSIK
                             ***************:****** ****:********:*.****** *****:*****.**

NODE1_KJHFFCBH_00001         RALSMVQEVPILGAETWYIDGGRRQGKAARAAYWTNTGRWQVMEIEGSNQKAEVQALLMA
MN630242_AMGANPBD_00002      RALSMIKDAPILGAETWYIDGGRKLGKAAKAAYWTDTGKWQVMELEGSNQKAEIQALLLA
U11820.1_EANIPDKN_00002      RALSMIQDAPITGAETWYIDGSRKQGKAARAAYWTDTGKWQIMEIEGSNQKAEVQALLLA
                             *****:::.** *********.*: ****:*****:**:**:**:********:****:*

NODE1_KJHFFCBH_00001         LKEGPEEMNIITDSQYILNIMNQQPDLMEGIWQEVLEEMEKKIAIFIDWVPGHKGIPGNE
MN630242_AMGANPBD_00002      LKAGPEEMNIITDSQYVINIILQQPDMMEGIWQEVLEELEKKTAIFIDWVPGHKGIPGNE
U11820.1_EANIPDKN_00002      LKAGSEEMNIITDSQYILNIINQQPDLMEGLWQEVLEEMEKKIAIFIDWVPGHKGIPGNE
                             ** * ***********::**: ****:***:*******:*** *****************

NODE1_KJHFFCBH_00001         EVDKLCQTMMIIEGDGILNKRPEDAGYDLLAAQEIHILPGEVRIVPTRTRIMLPKGYWGL
MN630242_AMGANPBD_00002      EVDKLCQTMMIIEGDGILDKRSEDAGYDLLAAKEMHLLPGEVKVIPTGVKIMLPKGYWGL
U11820.1_EANIPDKN_00002      EVDKLCQTMMIIEGEGILDKRSEDAGYDLLAAQETHFLPGEVRIVPTKTRIMLPKGHWGL
                             **************:***:** **********:* *:*****:::** .:******:***

NODE1_KJHFFCBH_00001         IMGKSSIGNKGMDVLGGVIDEGYRGEIGVIMINLSKKSTTILEKQKVAQLIILPCKHESL
MN630242_AMGANPBD_00002      IIGKSSIGSKGLDVLGGVIDEGYRGEIGVIMINLSRKSITLLEQQKIAQLIILPCKHEVL
U11820.1_EANIPDKN_00002      IMGKSSIGSKGMDVLGGVIDEGYRGELGVIMINLTKKSITILEKQKVAQLIILPCRHESL
                             *:******.**:**************:*******::** *:**:**:********:** *

NODE1_KJHFFCBH_00001         EQGEIIMNSERGEKGFGSTGVFSSWVDRIEEAELNHEKFHSDPQYLRTEFNIPRIVAEEI
MN630242_AMGANPBD_00002      EQGKVIMDSERGDKGYGSTGVFSSWVDRIEEAEINHEKFHSDPQYLRTEFNLPKMVAEEI
U11820.1_EANIPDKN_00002      QQGEIQMDSERGEKGFGSTGVFSSWVDRIEEAELNHEKFHSDPQYLRTEFNLPRIVAEEI
                             :**:: *:****:**:*****************:*****************:*::*****

NODE1_KJHFFCBH_00001         KRKCPLCRIRGEQVEGKLKIGPGIWQMDCTHFNGKIIIVAIHVESGLLWAQIIPQETADC
MN630242_AMGANPBD_00002      RRKCPVCRIRGEQVGGQLKIGPGIWQMDCTHFDGKIILVAIHVESGYIWAQIISQETADC
U11820.1_EANIPDKN_00002      KRKCPLCRIRGEQVGGQLKIGPGIWQMDCTHFNGKIIIVAVHVESGFLWAQIIPQETAEC
                             :****:******** *:***************:****:**:***** :***** ****:*

NODE1_KJHFFCBH_00001         TLKAIMQLVSTHNVTEIQTDNGPNFKNQKIEGLLSYMGIKHKLGIPGNPQSQALVENANS
MN630242_AMGANPBD_00002      TVKAVLQLLSAHNVTELQTDNGPNFKNQKMEGALNYMGVKHKFGIPGNPQSQALVENVNQ
U11820.1_EANIPDKN_00002      TVKALLQLICAHNVTELQTDNGPNFKNQKMEGLLNYMGIKHKLGIPGNPQSQALVENANN
                             *:**::**:.:*****:************:** *.***:***:**************.*.

NODE1_KJHFFCBH_00001         TLKVWIQKFLPETTSLDNALALALHCLNFKQRGRLGKMAPYELYTQQESLRIQDYFSQLP
MN630242_AMGANPBD_00002      TLKAWIQKFLPETTSLENALALAVHCLNFKQRGRIGGMAPYELLAQQESLRIQEFFSKIP
U11820.1_EANIPDKN_00002      TLKAWIQKFLPETTSLDNALALALHCLNFKQRGRLGKMAPYELYIQQESLRIQDYFSQIP
                             ***.************:******:**********:* ******  ********::**::*

NODE1_KJHFFCBH_00001         QKLMMQWVYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIRRVPEPCTLPE
MN630242_AMGANPBD_00002      QKLQAQWIYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIKRVPEPCALPE
U11820.1_EANIPDKN_00002      QKLMMQWVYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIRRVPEPCTLPE
                             ***  **:*****************************************:******:***

NODE1_KJHFFCBH_00001         GDE
MN630242_AMGANPBD_00002      GDE
U11820.1_EANIPDKN_00002      GDE
                             ***

Overall, the level of homology is quite high at the protein level. There's maybe a little evidence that this contig's protein is closer to U11820.1 than to MN630242 (though there are also several positions where U11820.1 is the outgroup). So I'm not sure why pieces of this contig are so hard to align at the DNA level without further investigation, but honestly overall of these are looking like very similar viral sequences, as might be expected.
I also aligned the big contig to both U11820.1 and MN630242 at the DNA level in 3-way clustalo alignment it is ok. Not sure what the difference is there. There is no obvious reason to prefer one reference over the other, according to my eye. When I have clustalo output a clustering solution the big contig is essentially equidistant between the two references, so I don't buy that U11820.1 is "better" as a reference:
(
MN630242.1:0.116151
,
(
NODE1_rc:0.0504856
,
U11820.1:0.0504856
):0.0656651
)
;

As I noted in the comments, it's unfortunate that QUAST and SPAdes aren't working very well. I don't know what the deal is there. But it seems like overall the assembly is quite close to both references. I don't have any intuition for why SPAdes is doing what it's doing, or why alignment at the DNA level is acting weird (lots of artificial gaps?). But I don't think it's a pathology of SPAdes particularly, it seems like it's just something weird that happened with this dataset interacting with DNA aligners. Possibly something went strangely with read deposition?
Update:
See OP's self-answer as well (and also discussion in the comments).

Based on the (extremely informative!) follow-up updates that you have provided, I think that we can tentatively answer your questions:
The tl;dr is that you can check your reads to make sure they aren't weird and you can check your output sequence to see if it's weird.

Am I using SPAdes correctly? Is there a better way to use it?

I think that you are using SPAdes correctly, based on your experience getting ~finished viral genomes in the past with the same workflow. (THis is also my experience using SPAdes on viral genomes with Illumina PE data, is that it "just works".)
It is possible that you are using some other tools in your workflow in a suboptimal way. For example, maybe your trimming is not quite correct (it does seem based on your answer that you are already trimming). The question in my mind is then whether there is some non-intuitive thing in the reads that should still be trimmed.
I would suggest running FASTQC on your reads for this run and also your reads for the other runs that have worked, to see if something weird is going on. FASTQC directly checks for adapter sequence, weird GC composition, quality scores etc. That will give you a lot of information.

Is it possible that the sequences that blast finds matching with the SPAdes contigs are indeed sequences that someone got using SPAdes and that's why these match (not because these are correct)?

This is formally possible, of course. Our main way of knowing that you're getting the right genome is... does it look like other genomes that have been generated in the same way?
If you are really worried about this, I would suggest directly inspecting the sequence to see what the matches are. If you like, you can then blast the matches to see if they
Another approach would be to annotate the genome and make sure that it "looks like" a phage genome. E.g. it has the expected genes from FIV. Prokka is very easy to run and works ok on phage.
A slightly self-serving suggestion is to use PhageTerm to check your reads against your final genomes to ensure that they behave like phage- e.g. they have termini in the correct places, etc. (it does this with test coverage.) I am a contributor to PhageTerm, so I happen to know it a bit. Likely there are other options.

What other option could explain these circumstances that I could be missing?

As I've suggested, I think that that leaves the reads. How sure are you that they are comparable to your other runs? They look like they should be sufficient from my glance at SRA, but maybe someone bumped the MiSeq while it was running. I've given suggestions above for checking with FASTQC.


maftools tcgaCompare plot: x-axis labels getting clipped

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


Does picard markduplicate toggle PCR duplicate samflag

Yes, if MarkDuplicates encounters a pair that's marked as a duplicate that it considers (for whatever reason) to not be a duplicate then it will unset the duplicate mark. You can test this yourself by making a small BAM file either with or without duplicate entries but some marked regardless. In the case of actual duplicates, the duplicate flag will be maintained. In the case where there are no actual duplicates then the flag will be removed. At least that's the case in version 2.14.


BAM to gene expression matrix (UMI counts per gene per cell),10X

You should be able to parse out what you need using the tags in the .bam.  10xGenomics' website says what tags they add.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam

Going backwards would also involve parsing the tags, because you have to make two fastq files, and a simple bam -> fastq pipeline won't do that correctly.


Why do BLASTn and prokka not seem to be searching the whole fasta file?

With a k-mer size of 28 it shouldn't be finding that many matches. And the prokka results are suspicious as well. Maybe you have multiple contigs (none larger than 100kb) in that file? What is the result of

grep ^'>' fasta_file | wc -l

? This would show how many contigs you have in the file.


How to make work programs from the $PATH?

echo $PATH will show the directories in your search path. Either copy the samtools executable there or add its location to the $PATH variable. You would be well served with a unix command line course or tutorial.


What is the standard way to work with a diploid reference genome? Complementary strands?

At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome.

The premise of the OP's question is false.  The natural reference representation of the human genome is not diploid.

Think of a reference genome as a map, and not as a specific example of a human being's DNA.

Not only is the human genome reference haploid, but it is also a composite genome. This means that the human genome reference sequence is composed of sequences from multiple individuals.  In other words, the human reference does not correspond to any one human sequence.

Any particular read from a DNA sequencer will be a read from a human genome that diverges from the reference genome. So an algorithm that tries to match the read to the reference genome will always need to handle potential discrepancies. Adding a second map against which to match a read would not change that fact.  Therefore, there is little value in providing a second haploid reference genome.

Side note:
There are parts of the human genome that "are too complex to be represented by a single path", and the Genome Reference Consortium provides "alternate loci" for such regions of the genome.

Answer to the OP's second question

  More importantly, how have large-scale genomic studies dealt with the fact that the haploid reference genome is a consensus-based "half" of a human genome?

I interpret this question as "How do large genomic studies represent genome diversity that cannot be represented by a haploid reference genome?"

One standard approach that has been very popular in genome-wide association studies is the use of a haplotype reference panel such as for example from the 1000 Genomes Project.  Modern phasing programs exploit the shared ancestry of the samples in a study and reference haplotypes to phase a study's samples. One popular phasing program is Impute2


How to remove duplicates from a fasta file using python

from Bio import SeqIO
import sys
import re

unique = []
fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile) > 250:
        gene_id = myfile.id
        if gene.id not in unique:
            unique.append(gene.id)
            mylist = re.match(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
            print (">"+mylist.group(1))

```


Are fgsea and Broad Institute GSEA equivalent?

According to the FGSEA preprint:

  We ran reference GSEA with default parameters. The permutation number
  was set to 1000, which means that for each input gene set 1000
  independent samples were generated. The run took 100 seconds and
  resulted in 79 gene sets with GSEA-adjusted FDR q-value of less than
  10−2. All significant gene sets were in a positive mode. First, to get
  a similar nominal p-values accuracy we ran FGSEA algorithm on 1000
  permutations. This took 2 seconds, but resulted in no significant hits
  due after multiple testing correction (with FRD ≤ 1%).

Thus, FGSEA and GSEA are not identical.

And again in the conclusion:

  Consequently, gene sets can be ranked more precisely in the results
  and, which is even more important, standard multiple testing
  correction methods can be applied instead of approximate ones as in
  [GSEA].

The author argues that FGSEA is more accurate, so it can't be equivalent.

If you are interested specifically in the enrichment score, that was addressed by the author in the preprint comments:

  Values of enrichment scores and normalized enrichment scores are the
  same for both broad version and fgsea.

So that part seems to be the same.


How to find the number of contigs produced? N50 length in base pairs? N90 length in base pairs?

With the data attached above, would my N50 and N90 length in BASE PAIRS so excluding the kb, be 7464 and 14314, respectively?

Looks like this is inverse so your N50 is 1032bp and your N90 is 106bp. N90 cannot be longer than N50 by definition. You L50 and L90 are 7464 and 14314

Would the number of contigs produced by this assembly be 7, since there are 7 different contigs produced or would it be 25,534?

25534
Now for the real answer. Do not calculate or report N50 on transcript assemblies. It is more than useless. Transcriptomes are - by nature - dynamic in length. Longer is not better in this regard. Are your short transcripts lower quality assemblies, probably not, but N50 treats them as such. Please read the Trinity wiki for other, better ways to QC your transcript assembly.


Index VCF file for rapid access by rsID?

I am not aware of any widespread index formats for query-by-QNAME for BAM or query-by-ID for VCF. I thought there were historical samtools{-devel,-help} threads to the effect of “in principle one could index a named-sorted BAM file, but there's insufficient demand so no-one's ever implemented it”, but I can't find any just now.

Moreover note that VCF files are by definition sorted by genomic position, so such an index would just be a hash table of rsID to file offset or genomic position.

If you are using dbSNP, probably you have a dbSNP database or API that you can use to query by rsID and get genomic position back. So you can retrieve your VCF contents by rsID in two stages: do that dbSNP query, and use its results to find the variants in your VCF file by position.


Infer the new (Ban et. al) ribosomal nomenlature (ex. uL53 ) from the ribosomal protein's sequence

This is a copy-paste nearly verbatim of comments so that the question has an answer

What you are trying to do is find what genes are in what homology cluster. This is common problem and there are many solution each with some issues.

Uniprot90 is indeed a cluster of homologues but it is too limited. Whereas you require clusters that span all the domains (universal).

COG

The nearly original grouping, is COG (cluster of orthologous genes). It is really handy for this but aren't really available from Uniprot and NCBI (its curator) has been phasing COGs out. In Uniprot records, they are patchily present as its derivative eggNOG (which has varying degrees of "scope", LUCA being the universal one):

<dbReference type="eggNOG" id="COG0451">
<property type="taxonomic scope" value="LUCA"/>
</dbReference>

Pfam

Another option that is well mapped and always present when possible in Uniprot entries is Pfam ids. These are domain family annotations and are also good option. It is a domain fold group. So domains within protein with homologous folds are grouped into one. If a protein has two known domains, it will have two Pfam entries —repeat proteins have loads, but most protein have one domain. In Uniprot XML you will find these are:

<dbReference type="Pfam" id="PF01370">

Similarly there is InterPro, which is more narrow than Pfam, but generally overlaps.

DIY

Then there's always the option of making your own psi-Blast motif on the NCBI website (setting high number of hits, setting the database to refseq or even PDB and doing a few iterations), saving the motif pattern (PSSM) and using local blast with it.

Caveats

In general you may however encounter these pitfalls:

you may have two protein clusters that share a single id
you have gene fusions (one protein, 2 domains of interest). In Pfam the protein will have two ids, with COGs or similar, you'd have only an anamalous length to go by.
you may have wished to find analogues (i.e. two protein that arose independently that have the same role (a lot of ribosomal protein are simply structural). No homology clustering scheme can find these.  Say a bacterium has evolved to have a different protein with that role, yet as it has never been studied it is a domain of unknown function (guilt by association is the common strategy, but that is beyond your scope here).
you have a very small protein which is poorly conserved


Plotting distance tree from blastn output

Specific questions

When there are multiple hsps in one alignment then which one do I include in the tree?

You can included them all because NJ is not computer intensive, so running loads is fine. You can remove identical sequences if you want either programmatically (these days I use pandas Python) or something like RAxML. Its honestly not critical for an NJ tree (and in fact the more the better) providing you can read the labels. A tree viewing program such as FigTree allows you to collapse clades to make the tree more readable to users.

If two sequences map to two different parts of the query sequence and do not overlap, then what should I do?

Hmmmm.... good question. Do 3 alignments one with one e.g. 5' fragment, the other with the 3' fragment and attempt a joint alignment. Under other non-pairwise tree buidling algorithms this would be disasterous, NJ gets around the problem because it is pairwise. HOWEVER, ask yourself whether your feel the resulting tree is meaningful, I would argue not because there are no homologous sites between the 5' and 3' partial sequences and phylogeny models point mutations. NJ will give you a tree however and will even bootstrap it.

Generic questions, you can actually do this as one of the new features of NCBI's Blast

Go to Blast here, https://blast.ncbi.nlm.nih.gov/Blast.cgi
Enter your sequence into the box (it doesn't accept PDB codes alone)
Enter the protein database - when I first did this calculation I used SwissProt, thinking there would be alot of sequences - I then used "nr"
Under the algorithm parameters enter "50" (default is too many)
Hit "Blast"
Once the search is complete at the top of the page are the hyperlinks: "Other reports:
Search Summary [Taxonomy reports] [Distance tree of results]"
Click on "Distance tree of results"
The tree page will load automatically, automatically aligning your sequences and producing, in this case a parsimony based tree, but there is also the option of a nj tree (recommended) ....
Click "Tool", "Download", "PDF" ...

You can generate a distance and matrix and/or perform a neighbor-joining phylogeny in MEGA10 and using Timura 3-parameter model is a good model. The bootstrap option is important. Its very easy to use, the format is a bit weird

# MEGA
Title: Stuff
# Sequence1
ACTAGACGT
# Sequence2
ACCTTAGGA

etc ...


Is it OK to use Blast+ to query NCBI's 16s rRNA database with 16s DNA sequences?

Your confusion is understandable, since it isn't explained in the NCBI description of the database, but most of the sequences in the database are derived from DNA that encodes the 16S rRNA gene, not from the rRNA itself. In some cases, sequences are derived from cDNA that is synthesized using the actual rRNA as a template, but it's still DNA sequencing. Direct RNA sequencing is extremely uncommon. Depending on the cDNA sequencing protocol (first strand vs. second strand sequencing), the transcript sequence might be the reverse-complement of the 16S rDNA, but the only difference that you would see is that the numbering of the bases would be reversed. If you are sequencing rDNA from organelles, there may also be a few bases that are different, due to RNA editing.


Print specific columns in a matrix on the basis of sample id's in the header

There's a million different ways to do simple taks like this. Pick one and learn it:) Python is my go to as it can do pretty much everything you'd ever want. I saved your example as csv with no spaces called expression.txt.
import pandas as pd
df = pd.read_csv('expression.txt')
df[['SM-5GZZ7','SM-5GIEN','SM-5EGGH']]

SM-5GZZ7    SM-5GIEN    SM-5EGGH
0   0   0.0 0.0
1   187 109 143 251 113 139 199 NaN NaN
2   0   0.0 0.0
3   1   0.0 0.0
4   0   0.0 0.0
5   0   0.0 1.0


where to find pathogenicity metadata for E.coli genomes

you can use BacDive. This is the largest database for bacterial metadata.:
https://bacdive.dsmz.de/
e.g. you can search for the sequence accession "NCTC8196" and you would get the metadata of the according strain:
https://bacdive.dsmz.de/search?search=NCTC8196

What kind of data is available, is highly dependent on the strain. Usually type strains are more rich on data that others.

For further questions don't hesitate to contact us:
contact@bacdive.de


I'm not familiar with the program, but apparently Hail is setting itself up as a swiss-army chainsaw project for doing downstream analysis on variant-called datasets.

An overview of Hail can be found here:

http://blog.cloudera.com/blog/2017/05/hail-scalable-genomics-analysis-with-spark/

A tutorial on association testing can be found here:

https://hail.is/hail/tutorial.html#Association-testing


how to calculate Pairwise Alignment Scores for blosum62

Your alignment is

SDRVIKAAIFDPIQPDF---G-----PVYFGLGHVH
RDLVERLFILDMI-PGLIKAGDSFPIPVALMINHIF

For each position in the alignment you calculate the score for that alignment. For position 1 we'd look up S vs R in the matrix and find a score of -1. You continue doing this until you hit the first -, which is not in the matrix. The first - is a gapopening, each consequent - in a series of -'s counts as a gap extension, instead of an opening.

---G--
IKAGDS

gives
-12 -2 -2 6 -12 -2


BLASTN should definitely report alignments in the order that the query sequences were provided.

Based on your description and my personal experience, I think the most likely explanation is some kind of mix up. It happens to all of us, especially when we're in a rush to get an answer, but even sometimes when we're being careful and disciplined.

If you can go back and recreate the result, I would be very interested to see the two pairs of input and output files. If not, I wouldn't spend too much more time worrying about it. :-)

UPDATE Alignments are reported in the order that the query sequences are provided, even in tabular mode. The order of each query's hits is sorted by bitscore or e-value, but that is a secondary sorting. I've included an example below that has two queries, seq1 and seq2 with -outfmt 6.

seq1    scaffold2       100.000 280     0       0       1       280     11131   11410   7.73e-146       518
seq1    scaffold294     88.693  283     25      4       1       279     12160   12439   6.46e-92        339
seq1    scaffold1525    85.315  286     29      9       1       277     27142   26861   3.07e-75        283
seq1    scaffold1874    83.505  291     32      10      1       280     8558    8273    1.86e-67        257
seq1    scaffold147     82.653  294     34      9       1       279     10797   10506   1.45e-63        244
seq1    scaffold2478    82.374  278     38      8       11      280     7622    7896    1.13e-59        231
seq1    scaffold3405    87.940  199     17      5       1       196     14390   14584   1.46e-58        228
seq1    scaffold395     83.060  183     23      4       106     280     12021   11839   5.40e-38        159
seq1    scaffold2853    85.965  114     13      3       1       113     20864   20753   9.16e-26        119
seq2    scaffold9       100.000 140     0       0       1       140     1121    1260    2.34e-68        259
seq2    scaffold2950    90.909  110     6       3       1       106     110     1       6.84e-34        145
seq2    scaffold3416    95.385  65      3       0       1       65      11882   11818   1.16e-21        104
seq2    scaffold3103    95.385  65      3       0       1       65      9736    9672    1.16e-21        104
seq2    scaffold5297    88.608  79      7       1       64      140     150     228     6.99e-19        95.3


Is there a tool that can perform a read-group-aware mpileup from a single file?

I broke down and finally scratched this itch. I have implemented a tool that performs read-group-aware text-pileup from a single SAM/BAM file. The tool is called streaming_pileupy and it's available from pypi for installation with pip. After installation the command would be:

spileup input.bam sampel_names.txt

The tool is pre-alpha and missing some features. Contributions are welcome.