Why is bwa-mem the standard algorithm when using bwa?

bwa mem is newer, faster, and [should be] more accurate, particularly for longer reads.

From the bwa man page (presumably in Heng Li's own words):

  BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. It consists of three algorithms:  BWA-backtrack, BWA-SW and BWA-MEM. The first algorithm  is designed for Illumina sequence reads up to 100bp, while the rest two for longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar  features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and  more  accurate. BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.


Run Nextflow with file dependencies inside a Docker container?

The following workaround seems to work but ideally I don’t want to have to change the script portion of my process so that it works with Docker.

params.ref = 'data/reference/GRCh38.fa'

ref = file(params.ref)
ref_index_dir = ref.parent
fastq_channel = Channel.fromFilePairs('data/fastq/*_R{1,2}.fastq.gz')

process map_reads {
    input:
    set sample, '?.fastq.gz' from fastq_channel
    file ref
    file ref_index_dir

    output:
    file("${sample}.sam") into sam_channel

    script:
    """
    echo 'bwa mem -Y "reference/$ref" 1.fastq.gz 2.fastq.gz >${sample}.sam'
    stat -t reference/${ref}.bwt && touch ${sample}.sam
    """
}

Here, reference inside the script section refers to the parent directory name of the reference, and needs to be adjusted accordingly. Having to refer to reference/$ref is annoying but at least it works — both inside a Docker container and outside of it.


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

import pandas as pd

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

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

etc.

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


MultiQC can merge all your different reports into a single one.
Which could be useful once you manage to know which QC tools to use.


How can I identify a recessive and dominant gene?

I'd probably look up the genes in OMIM, which provides detailed description of the functions and phenotypes associated with genes, including inheritance mode.

For example, this is what you would get for TP53:


Fitting a principal curve into a diffusion map

If fit$ord gives you the order of cells on the curve, you should be able to obtain the first and last cells by doing the following:

first_cell_index <- which.min(fit$ord)
last_cell_index <- which.max(fit$ord)

However, I'm not sure this is so informative - if you're looking to get representative gene expression for the start and end of the curve, you should really be looking at the cluster of cells towards the start and end of the curve - for example, by comparing cells in clusters 1 and 3 that you showed in your first plot. A single cell will suffer from dropout and not tell you a lot about the cluster.


Error in formation of correlation matrix-(corrupt matrix -- dims not not match length)

You can't compute a correlation matrix when you have more than ~46000 rows, since a standard R matrix can have a maximum of 2^31-1 values. Have a look at packages like bigcor. Alternatively, consider if having 3+ billion correlation values is really useful.


ANY (technical) reason behind submitting sequences to GenBank versus ENA Sequence

Opinion: I prefer to use EBI/ENA because the turnaround time for "customer service" on editing metadata about a submission is faster than on NCBI. Inquiries to change metadata on NCBI can go unanswered for up to six months in my experience.

Answer to your question: No, there are no technical reasons that make it preferable to submit data through one website or the other. All entries end up being hosted by both databases anyway. The preference is up to the user to select whichever upload interface they prefer.


Obtaining identifier from plot of sequence GC%

The tricky bit here is that you need a list of sequence records sorted in the same order as the CG values. You can accomplish this by making tuples of the seq records and GC values, then sorting the tuples specifying the GC values as the sorting key:

gc_w_records = sorted(((GC(rec.seq), rec) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")), key=lambda tup: tup[0])
gc_values = [x[0] for x in gc_w_records]
records = [x[1] for x in gc_w_records]

I don't know much about matplotlib and it's going to be challenging to display 100 IDs, but you'll probably be interested in the extremes of the GC distribution, which you can obtain directly from the sorted list of records:

records[0:2]
# [SeqRecord(seq=Seq('CGTCACGAGGTCTCCGGATGTGACCCTGCGGAAGGATCATTGTTGAGATCACAT...CAT', SingleLetterAlphabet()), id='gi|2765587|emb|Z78462.1|PSZ78462', name='gi|2765587|emb|Z78462.1|PSZ78462', description='gi|2765587|emb|Z78462.1|PSZ78462 P.sukhakulii 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...AGG', SingleLetterAlphabet()), id='gi|2765568|emb|Z78443.1|PLZ78443', name='gi|2765568|emb|Z78443.1|PLZ78443', description='gi|2765568|emb|Z78443.1|PLZ78443 P.lawrenceanum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]

or:

records[-3:-1]
# [SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAGCAG...TAA', SingleLetterAlphabet()), id='gi|2765656|emb|Z78531.1|CFZ78531', name='gi|2765656|emb|Z78531.1|CFZ78531', description='gi|2765656|emb|Z78531.1|CFZ78531 C.fasciculatum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGCTCATTGTTGAGACCGCAA...AAG', SingleLetterAlphabet()), id='gi|2765625|emb|Z78500.1|PWZ78500', name='gi|2765625|emb|Z78500.1|PWZ78500', description='gi|2765625|emb|Z78500.1|PWZ78500 P.warszewiczianum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]


Can I use a regular liner regression model when I'm working with DNA methylation data?

You can use either, but lmFit has the benefit of returning an object that can be used with eBayes() so you can pool information across genes/probes/whatever. lm() is a base R function applicable basically everywhere. lmFit() if from the limma package, so originally intended for microarray data, though these days pretty much everything omics is analyzed with it since sample sizes tend to be small so everyone want's to pool information for more statistical power.

BTW, if you have methylation data then you'll want to fit data in logit space.


Bioconductor packages are R packages that can only be used directly from python via rpy2 or similar. If you've ever had reason to use rpy2 you'd probably realize quickly that while it works, it's no where near as easy as doing the same thing directly in R.

Biopython is a package for doing a number of useful things, but it can in no way come remotely close to performing even a percent of what Bioconductor packages do. That would be an impossible feat for a single python package.


Any suggestions for cultivar identification using SSR (simple sequence repeat) markers

Use Fst (migration statistic) such as the Fstat program to generate a distance matrix and solve it via clustering, such as UPGMA or neighbor-joining. This analysis is particularly useful for diploids. I presume you have heterozygotes in your analysis. You can switch Fis (inbreeding) if you have excessive homozygosity.
The other way in is a haplotype network. When I get a minute I'll post the link to the programs. Machine learning isn't needed in this instance, although you could do a PCA style analysis, i.e. unsupervised learning. Unsupervised learning doesn't need a large sample size.

The other solution is a stepwise mutation phylogeny, but this is for microsatellites (one to four repeats). This is called Delta Mu squared and describes the method of repeat unit mutation. It is a distance method.
Again when I have a moment I will supply weblinks.


What are the standard ways to visualize protein-protein or gene-gene interactions?

since I am not yet allowed to comment, I will have to pots this as an answer.

I guess Cytoscape (http://cytoscape.org/) would be a nice way to analyse such interaction lists/networks (still depends a bit on what you plan to do).

However, for this you will have to preprocess the downloadable excel sheet a bit.

save as tab-separated file (tsv)
run the below python3 script to format into an easily importable format for cytoscape. You can run the script as python3 script.py <path-to-tsv-file> <delimeter for cytoscape> . The new format is written to std out on the console, therefore maybe you need to write the output to a new file. A sample call could look like python3 testbench.py human_interactome_2017-06-12.tsv.csv ";". Make sure to add the delimeter in quotes.
you can then import the new file into cytoscape (import network as table) where the first column is your source, the third column your target and the middle column your interaction type.
Within cytoscape you can then adapt your edge style to vary the line type in a discrete mapping depending on the value of the interaction.

Screenshot from cytoscape:

Code below:

import sys

HIIIidx = 4
HII5idx = 5
venkatesanidx = 6
yuidx = 7
Litidx = 8

print(sys.argv)

if len(sys.argv) <= 1:
    print("Usage: script.py file delimeter")

delim = ';'
inFile = sys.argv[1]

if len(sys.argv) > 2:
    delim = sys.argv[2]

def printInteraction(syma, symb, inter):
    print(syma + delim + inter + delim + symb)

with open( inFile , 'r') as file:

    file.readline()

    for line in file:

        aline = line.strip('\n').split('\t')

        if aline[HIIIidx] == '1':
            printInteraction(aline[1], aline[3], 'HIII')

        if aline[HII5idx] == '1':
            printInteraction(aline[1], aline[3], 'HII5')

        if aline[venkatesanidx] == '1':
            printInteraction(aline[1], aline[3], 'VEN')

        if aline[yuidx] == '1':
            printInteraction(aline[1], aline[3], 'YU')

        if aline[Litidx] == '1':
            printInteraction(aline[1], aline[3], 'LIT')


You made the hardest work, you were not far to achieve what you wanted ;)
This script.py should work after you made it executable (in a terminal type: chmod u+x script.py), just type ./script.py in your terminal in your working directory (so here you should be in your "sp" directory). Note that I got inspired by: https://www.biostars.org/p/42474/ to know how to launch pymol within the script and also that I would personally not use pymol for doing this kind of analysis/data extraction from a pdb/pdbqt file but rather use a dedicated analysis package (there are plenty like MDanalysis). I personally use pymol (pymol script) for visualization purpose.
#!/usr/bin/env python

import __main__
__main__.pymol_argv = ['pymol', '-qc']  # Quiet and no GUI

from pymol import stored
import pymol, os

pymol.finish_launching()

#put what you've done into a function
def myfunc(file):
    spath = os.path.abspath(file)
    sname = spath.split('/')[-1].split('.')[0]

    pymol.cmd.load(spath, sname)

    stored.a, stored.b, stored.c, stored.d = [], [], [], []
    pymol.cmd.iterate_state(1, "%s///UNL`1/N*" % sname, "stored.a.append(ID)")
    pymol.cmd.iterate_state(1, "%s///UNL`1/O*" % sname, "stored.b.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASP`113/OD2" % sname, "stored.c.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASN`312/ND2" % sname, "stored.d.append(ID)")

    f = open('%s.txt' % sname, 'w')

    for a in stored.a:
        for c in stored.c:
            d1 = pymol.cmd.distance("id %s" % a, "id %s" % c, cutoff=5)
            f.write("%s D113 %s\n" % (a, d1))
            print(d1)

    for b in stored.b:
        for d in stored.d:
            d2 = pymol.cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
            f.write("%s N312 %s\n" % (b, d2))
            print(d2)

    f.close()

mypdbqt = os.popen('ls *.pdbqt').read().split()

for i in mypdbqt:
    myfunc(i)

pymol.cmd.quit()


Agreed with the above I would use an existing dataset from NCBIs dbSNP here
Random doesn't work well because HW signal will be meaningless as soon as the Jukes-Cantor distance is exceeded (which it will) and plink needs homologues to generate any signal either in 'reasonable' panmixia or clonality. The result will be panmictic beyond even the most panmictic organism ever known. I don't mind having a simulation of 100% clonalty, because 100% clonality does occur, either through inbreeding or binary division, but this simulation would generate the complete opposite result.


Putting cells in Seurat dimplot selectively

If you're using a GUI you could select the cells interactively:

plot <- DimPlot(seurat_obj, reduction = "umap")

Then select the cells by clicking around them

select.cells <- CellSelector(plot = plot)
Idents(seurat_obj, cells = select.cells) <- "SubCells"

and subset based on these cells

sub_cells <- WhichCells(seurat_obj, idents = "SubCells")
sub_obj <- subset(seurat_obj, cells = sub_cells)

You could then use subobject in your downstream dimplot.

See here for some of the interactive options. I guess you could outright subset on UMAP coordinates, but that may not always be easy if your boundary isn't linear.


Can someone recommend software for reference-based genotype-estimation?

In case anyone was interested, a new piece of software for rapidly imputing and phasing low coverage sequence data from genotype likelihoods has just been released. It is called GLIMSPE and is exactly what I was after.

It looks to be much faster (1200x speedup) compared to Beagle 4.0, and is more accurate at imputing low frequency variants. It is also more amenable to use with larger reference panels (such as the Haplotype Reference Consortium panel).


Recommendations for missing value imputation - DNA methylation data

I think this is still an active field of research.  I have heard of Phenix, which might be appropriate.


GATK best practices are explicably meant to consume BWA MEM generated BAM.  Whilst SNAP may be faster, the Broad will not have tested it for compatibility with GATK as such you can't guaranty using it won't have unexpected consequences.

As such you'd be better off using BWA MEM because I assume accurately called variation is always better than fast and incorrectly called variation.  The main issue you'll have is ensuring shorter split hits and mapping quality are reported in the same way as bwa MEM -M which GATK/Picard is expecting.  Ultimately however you'd be better off posting this question on the GATK forum.

It's also worth noting that the soon to be released GATK 4 will utilise bwaspark which can distribute it's alignment processes across Apache Spark for increase performance.  Consequently I can't see SNAP being adopted anytime soon.


How can I not show insertions in the Integrative Genome Viewer (IGV)?

Preferences -> Alignments -> Hide indels < 'x' bases.

Credit to Pierre Lindenbaum


What is the best method to estimate a phylogenetic tree from a large dataset of >1000 loci and >100 species

This paper claims that FastTree is almost as accurate as RAxML, while being much faster. You just have to be careful, however, that the support values output by FastTree are not bootstrap values, they are based on the Shimodaira-Hasegawa test. (Also, see this comment for the case you have very short branch lengths).   [update: However, according to the recent comparison paper mentioned below FastTree performed quite poorly in comparison to RAxML or IQ-tree.]

From what I understand, you should use ExaML only if your data is too large to be handled by RAxML in a single node. ExaML should perform like RAxML but with some parallelization overhead. For all effects I treat them as the same. I don't know of relevant advantages of phyML over RAxML (for me, it's easier to use but I am very used to phyML).

I am not familiar with IQ-tree, but its authors claim that even given the same time as RAxML or phyML, IQ-tree already finds better likelihoods more often than not (although by default it takes a bit longer to converge). A recent comparison between all these programs favoured IQ-TREE for both single-gene and concatenation analyses (with RAxML very close). It may also estimate branch support through a SH-like test only, but I'm not sure.   [update: IQ-tree offers 3 measures of support, standard bootstrap, aLRT, and ultrafast bootstrap. See OP's comment below for details.]

However, since you have few missing data, you might also want to try a single-locus tree inference followed by gene tree clustering (using treescape or treeCL) to see how spread your data is, or to see the effect of removal of outliers, or to use ideas similar to statistical binning.


Which one is a more convenient assembly?

I asked that on Reddit. I didn't get an answer to all my questions but I got an answer for how to assemble the overlapping contigs.
If I use Flye I can assemble these and get larger contigs, better NG50, and higher genome fraction.

s-aligner + Flye assembly of SRR8357422 Genome fraction:95,95% - Largest Alignment: 134.103 - NGA50: 134.103 - LGA50: 1
Given that the data in the first assembly apart from being more useful itself in my opinion, is objectively more suitable to generate a better result in all senses with a simple additional step that takes a few seconds, my only possible conclusion is that it is overall a better result and therefore this sort of result will often be preferable.
Of course, there are other aspects to consider, like error rate, that are not considered here as it wasn't part of the question (mainly because I think Quast is not adequate for calculating it). But BTW, mismatches per 100kbp are also lower.
[EDIT] I want to add that this is not granted that in all cases Flye will generate a clearly better assembly.


Looking for a way to download all pdb files that detail a protein-drug complex

Proving a negative is hard, but to my knowledge I do not know a database that has all the structure of drug bound protein. I can give some close matches.

First a few terms to be on the same page. A ligand is something that binds, which can be physiological, such as a cofactor or substrate, or not (especially man-made), such as a drug, i.e. a compound that has a physiological effect, or a hit, a small fragment from a fragment screen that will not bind with nanomolar affinity. Generally, ligands are small molecules, but may include oligopeptides and oligonucleotides. In the case of crystal structures, you also have crystallisation additives, which are not ligands, but are hetatoms in PDBs. Also, the definition of drug may vary. drugbank.ca has the following categories of drugs "approved", "experimental", "withdrawn" for pharmaceuticals, "illicit" for narcotics and "investigational" for probes and counts macromolecules as drugs.

RCSB PDB: drug mapping table is a simple searchable table with ligands that are clinical drugs only —experimental compounds (=probes) will be absent.
bindingdb.org contains many table, including a '3D' one. This table contains all the PDBs with ligands and their affinities. This contains crystallisation additives and natural ligands. Peter Curran's repo for the HotSpots algorithm (first written by Chris Radoux) has a blacklist, but you are bound to get stuff missed.
bindingmoad.org is a database with all the protein-ligand structures from PDB up until 2018, but it includes natural ligands, like cofactors and substrates as above.
https://klifs.vu-compmedchem.nl a well curated database of kinases, which has a lot of members and many ligands, but is not at all diverse. As a result it is a clean dataset often used as a testbed for novel methodologies.
I mentioned in the comments the SGC Target Enabling Packages as a small subset of proteins with many hits each. But I realise that there is not a handy PDB list and the targets are not drugs, but probes.


The CART tool let's you upload a set of names and map them (optionally in a fuzzy way) to STITCH 4 identifiers, and then use those to map to ATC codes (using the chemicals sources download file). It's a bit indirect, and I'm not sure what CART will do with the dosage info you mention.


How to remove all BAM read groups from all reads (not just the header)?

You might have to manually strip those auxiliary tags off:

samtools view -h your.bam | grep -v "^@RG" | sed "s/\tRG:Z:[^\t]*//" | samtools view -bo your_fixed.bam -

The sed bit is searching for the aux tag and removing everything up to the next tab.


The HIVdb page states the following in the second paragraph (emphasis mine):

  Sequences can be entered as plain text if just one sequence is entered. Sequences must be entered using the FASTA format if multiple sequences are entered. Sequences can be pasted in the text box or uploaded using the File Upload option. The upper limit is currently 1000 sequences containing 3000 nucleotides per sequence.

The sequence in your question is a protein, so it consists of amino acids, not nucleotides. So, if you instead give a nucleotide sequence, for example this one:

>KX692212.1 HIV-1 isolate CH47 from Switzerland pol protein (pol) and gag protein (gag) genes, partial cds
CCTCAGATCACTCTTTGGCAACGACCCCTCGTCACAATAAAGATAGGGGGGCAACTAAAGGAAGCTCTAT
TAGATACAGGAGCAGATGATACAGTATTAGAAGAAATGAGTTTGCCAGGAAGATGGAAACCAAAAATGAT
AGGGGGAATTGGAGGTTTTATCAAAGTAAGACAGTATGATCAGATACTCATAGAAATCTGTGGACATAAA
GCTATAGGTACAGTATTAGTAGGACCTACACCTGTCAACATAATTGGAAGAAATCTGTTGACTCAGATTG
GTTGCACTTTAAATTTT

The program will work as expected:

For future reference, I found the KX692212 sequence above by taking your protein sequence and running a tBLASTn analysis against HIV sequences in NCBI's nr database. That was the first hit.


Which tools are available to call somatic mutations on non-coding regions from whole genome sequencing data?

Those will have to be treated as separate analysis.  You can use Mutect or VarScan for the matched samples and the tools listed in the other question for the unmatched samples however that means that you won't be able to compare the results easily.  If the inter-sample comparison is important you will have to run all samples as unmatched.  Personally I have used LoFreq for somatic variant detection with some success.  It is very sensitive so the false positive rate can be high depending on how noisy your data is.


How to determine a correct promotor sequence? (Bisulfite sequencing - uni assingment)

I think there is no fixed size. Ideally you would have something like histone ChIP-seq or ATAC-seq data to approximate the relevant upstream interval. Since this is not always possible most people go with -500bp. Others might have different suggestions though. Maybe matching known motifs of transcription factors that regulate that promoter to the upstream region can help approximating a better window. But as Sanger sequencing usually gives about 1kb of sequence I'd just take a 500bp to 1kb window and sequence that, fwd and rev primers for sequencing to have double confirmation.


Seeking explanation of the hg38 files downloaded from bowtie 2 website

Bowtie is providing you with the index files for hg38.  This are the result of the bowtie2-build indexer.  However you have downloaded the indexes for Bowtie1 not bowtie2

  Bowtie 2’s .bt2 index format is different from Bowtie 1’s .ebwt
  format, and they are not compatible with each other.

Then according to the manual you pass the basename of the index for the reference genome. So this should be the name prior to the extension, .1.bt2 for Bowtie2.

bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --interleaved <i>} -S [<sam>]
Main arguments
-x <bt2-idx>
The basename of the index for the reference genome. The basename is the name of any of the index files up to but not including the final .1.bt2 / .rev.1.bt2 / etc. bowtie2 looks for the specified index first in the current directory, then in the directory specified in the BOWTIE2_INDEXES environment variable.


COVID-19 GWAS: what is a cross-replicating association?

Replication in GWAS is an analytical term, approximating to 'robustness'. It is not biological phenomenon of a genome replicating. What the authors are saying is these genes associate with an allelic variant(s) which correlate with the disease phenotypic of COVID-19.
Three of those genes are very interesting indeed notably CCR9, CXCR6 and XCR1, because these are chemokine receptors. The later for example is associated with regulating between the innate and adaptive (CD8+, i.e. cytotoxic T-cell) immune responses (described in the image below). The innate/adaptive interface is an extremely trendy area in immunology.
The authors are implying the immunological signalling associated with one or more of these chemokines is fundamental to combating COVID-19 infection. 'Replication' is simply the GWA geek for 'its real'.


If helpful for anyone else, I found that all and much more are available. See below 2 links.

https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/dataset.cgi?study_id=phs000424.v8.p2&phv=169091&phd=3910&pha=&pht=2742&phvf=&phdf=&phaf=&phtf=&dssp=1&consent=&temp=1

https://ftp.ncbi.nlm.nih.gov/dbgap/studies/phs000424/phs000424.v8.p2/pheno_variable_summaries/phs000424.v8.pht002742.v8.p2.GTEx_Subject_Phenotypes.var_report.xml


Are crystal structure of Ligands used in volunteer computational docking like FightAIDS@home?

I don't know that project's details, so cannot say officially, but undoubtedly it will be

crystal structures of protein or models of point mutants based on structures of protein.
Computer structure of the ligands based on force-field calculation

Protein structures

Firstly, HIV virus has been extensively studied by structural biologists, so most components have structures.

Secondly, docking on models is generally done only for threaded models based on closely related templates, say a human model of a mouse structure or a recently diverged pair (i.e. conserved proteins that split with the genome double duplication in bony fish). Whereas a given virus protein will likely have poor homology to a known structure and not be a good candidate.

A corner case is core point mutants... If you have a mutation in a protein it is likely going to be on the surface and not core: a model of a point mutant is likely very faithful to an empirical version —and generally you'd want to design drugs to a region that is not variable anyway as you'd want it to work universally. If you have a change in the active site to confer resistance say (i.e. in the core), there is generally a fitness loss to the organism (say lowered catalysis or thermostability) which needs to be compensated with other mutations (epistasis). As this is rarer and very interesting scientifically and medically, these resistant variants will most likely have been crystallised as the conditions for crystallization will most likely be the same.

Small molecule structures

Force field calculations for small molecules are most often valid (there computer Vs empirical RMSD is minute). The partial charges may be off, but that's not something crystals give anyway: In fact a ligand used in docking is (generally/ought to be) energy minimised beforehand. For docking getting charges wrong is the killer.

Crystallization of small molecules is used as validation of a synthesis route, so some may have structures, but they will not be used in a pipeline as a SMILES string is all you need.


CellPhoneDB is a publicly available repository of curated receptors, ligands and their interactions.

The interactions can be searched or downloaded.


How I can change the name of multiple files at once in R or terminal?

Use perl-rename. This is usually called rename on Debian-based systems like Ubuntu or Mint, and perl-rename or prename on others. Assuming you have it as rename, you can simply do:

rename -n 's/.*\.vcf/"A" . ++$c . ".vcf"/e' *snp.pass.vcf

The -n causes rename to only print what it would do, without doing anything. So once you make sure it does what you want, run the command again without the -n.


estimate genome size: kmer-based approach from PacBio reads

I don't think there is a method that would estimate a genome size using raw long reads.

The genome size estimates based on raw reads are done by fitting a model to kmer spectra (for instance Genomescope). The kmer spectra built from long reads are really messy due to the high error rate of long reads. That makes fitting of a model quite difficult. These methods assume small error rates. If the error rate is as high as 0.1 (which is normal for long reads), then the probability of having one kmer of 23 bases correct is ~8%. Meaning that if the coverage is 100x, just ~8x of homozygous kmers would be the real genomic kmers and it's very likely that kmer peaks of the sequencing errors, heterozygous kmers and homozygous kmers will just blend into one.

I can imagine, that if your coverage would crazy high (like 600x or more) the kmer spectra could separate errors and genomic kmers. You could also try to decrease kmer size to increase the proportion of correct kmers.

Another option would be to use a trick that they used to generate solid kmers for genome assembly using generalized de brujin graph using noisy reads. Maybe it would be better to take a look on the newer version of the assembler, that is called Flye.

But to be honest, I think for a long reads, a genome assembly is much more reliable way how to predict a genome size than any of the kmer tricks. I would simply assemble the genome, remap back reads, check the uniformity of coverage and conclude the genome size out it.

-- edit --

I tried to find some comprehensive evaluation of the assembly sizes compared to genome sizes, but I could not, so I just checked two examples, one big, one small and they both match quite well. Ammopiptanthus nanus was assembled to 823.74 Mb while kmer estimate was 889 Mb (ref). Plasmodium falciparum was assembled to 23.5 Mb while the original genome size estimate is 22.8 Mb. I know it's just two, but there is also a good reason to think that the genome size will be reasonable, feel free to add here more examples.

If you would be worried about contamination in the assembly, you can first run it though one of the methods to detect a contamination (like Blobology).


Does module size decrease as the number of samples increase in WGCNA?

The biggest difference a large number of samples makes is that you can usually decrease the soft thresholding power. If you use power 6 (or 12 for a signed network), try decreasing it to 3 or even 2 (6 or 4 for signed network). Apart from that, it could also be that the co-expression modules in your data are simply smaller than what people usually find, nothing apriori wrong with that.


How do kmer counters determine which kmer is 'canonical'?

When a k-mer is identical to its reverse complement, both are canonical. Note that a canonical k-mer is a sequence, irrelevant to its position(s) in the input string. More precisely, give a string $s$, its canonical string is
$$
{\rm canonical}(s|h)=\left\{\begin{array}{ll}
s & \mbox{if $h(s)<h(\overline{s})$}\\
\overline{s} & \mbox{otherwise}\\
\end{array}\right.
$$
where $\overline{s}$ is the Watson-Crick reverse complement of $s$ and $h$ is an arbitrary string hash function. In practice, we most often use 2-bit encoding as $h$. Order under such a hash function is equivalent to the lexicographical order.


Dmpfold error "No such file or directory: 'best_qdope.pdb'" when running "run_dmpfold.sh"?

The target variable, which is set to the first part of the sequence file name, should be sequence according to your given command line. However the error message thinks it is c600m2-tail-before-lysin-sequence. Have you edited the command line for simplicity?

I'm a little worried about the floating point errors, which are almost certainly from CNS as there are 3 CNS calls and the error shows up 3 times. I would test your CNS installation to see if everything is okay - maybe try some examples from their website. An error there would explain why no output is generated.

You can ignore any errors below line 5 of the above program output for now, as they are follow-on errors. In particular the question title isn't the 'real' error here, the real error is that no PDB files are generated on line 160.

By the way, I wouldn't use 10 500 as parameters for iteration number and number of models. It will take ages to run. I'd just use the defaults (omit the arguments), which we find gives good results.

Finally, I'll point out that the DMPfold repository has an issue tracker where you can raise any future issues. Though we may as well finish debugging this one here since we've started!


CCP4 file to a Python 3 numpy array or similar workaround

You can use Gemmi.

import gemmi
ccp4_map = gemmi.read_ccp4_map('my.map')
ccp4_map.setup()  # optional
np_array = numpy.array(ccp4_map.grid, copy=False)
# ...
ccp4_map.write_ccp4_map('out.ccp4')

The setup() call expands data to the whole unit cell.
The numpy array above accesses the data through Python's Buffer Protocol (it does not copy the data). So changing values in numpy array changes gemmi.Ccp4Map.

You can also use gemmi to read mtz file and transform map coefficients to a real-space map, see https://gemmi.readthedocs.io/en/latest/hkl.html.


Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose.

C++, Rust, Go, Nim and Julia can be as fast.

// Download https://raw.githubusercontent.com/lh3/minimap2/master/kseq.h
// Compile with: gcc -O2 -o myprog this-file.c -lz
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
// STEP 1: declare the type of file handler and the read() function
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r"); // STEP 2: open the file handler
    seq = kseq_init(fp); // STEP 3: initialize seq
    while (kseq_read(seq) >= 0) { // STEP 4: read sequence
        int i, sum, low;
        if (seq->qual.l == 0) continue; // no quality
        if (seq->seq.l < 15 || seq->seq.l > 30) continue;
        for (i = sum = low = 0; i < seq->seq.l; ++i) {
            int qual = seq->qual.s[i] - 33; // assuming Sanger quality
            sum += qual;
            if (qual < 20) ++low;
        }
        if (sum < seq->seq.l * 30.0 || low > 2) continue;
        printf("@%s\n", seq->name.s); // output can be made faster
        puts(seq->seq.s);
        printf("+\n");
        puts(seq->qual.s);
    }
    kseq_destroy(seq); // STEP 5: destroy seq
    gzclose(fp); // STEP 6: close the file handler
    return 0;
}

To apply the program to all *.fq in a directory (again, not tested):

mkdir -p filtered
ls *.fq | xargs -i echo ./myprog {} \> filtered/{} | sh

You can also replace the final sh with parallel, but the above should be fast enough.


TCGA bam file duplicates (aligned reads) and related last part of TCGA barcode

If you check the bam list of TCGA, you'll see that they use different numbers for the bam files (maybe different settings). What the difference is, I don't know, but if you want the file that they used for their paper, it is C500.TCGA-BL-A3JM-01A-12D-A21A-08.2.bam, the other one is not in the list.


Sequencing center returns data in several files

Usually, a single sequencing library is created from each sample. Assuming you are using barcoded indexes to identify multiple samples, the index is attached to each library, and the libraries are pooled to make a library pool.  A library pool can be frozen, and is very stable.

A pool can be sequenced multiple times by pipetting a portion of the pool and supplying it to the sequencing process.  Each time, a sample of the molecules in the pool will be sequenced.

Generally, it is true that sequencing runs are very repeatable, and that you will not suffer any batch effects.

If you want to check for batch effects you would need to analyse all of your sequence files separately, and check this first, or in parallel, to your main analysis.

One thing you might want to check is whether the same library pool was sequenced multiple times, or whether separate pools or libraries were made for each sequencing run.  In the second case, there may be more likely to be batch effects.


How to normalise scRNASeq data for differential expression analysis

I would suggest using a likelihood ratio test for differential expression using logistic regression with batch as a latent variable. In Seurat you can do:

markers <- FindAllMarkers(object, test.use = 'LR', latent.vars = 'batch')

(change "object" and "batch" accordingly)

See https://www.biorxiv.org/content/early/2018/02/14/258566 from Lior Pachter's group for some more info on the LR test.


Is there a way to measure cell line similarity using python?

At its current form, your question is not clear, for example which aspects of these cell lines are your trying to compare?

Regarding your point on the ontologies, in Python you can use the pronto library to work with ontologies in OBO format.


A site is simply an individual discrete position, normally a single nucleotide (or amino acid). For example, consider a toy alignment:

sequence1 ATGC
sequence2 ATCC
#           ^ arrow points to the third site

Here there is 1 mutation in 4 sites, so 0.25 substitutions per nucleotide site.

I'm not sure where the term originated, but probably one of the earlier definitions is from renowned population geneticist Motoo Kimura, who simply stated in 1968:

  In this paper, “site” refers to a single nucleotide pair

A pair in this context, means either a homozygous or heterozygous site.

Synonymous sites are positions where nucleotide changes lead to no amino acid change, these are typically at the 3rd position of the codon where the nucleotide is often degenerate, for example, TAT, TAC TAG, and TAA all encode the Serine amino acid. By contrast non-synonymous sites are where a single nucleotide change does lead to an amino acid change, and therefore possibly a more important functional change to the resulting protein.

Thus, the ratio of Ka/Ks (or Dn/Ds) is thus often used as a crude proxy of selection, with Ka/Ks < 0, Ka/Ks = 0 and Ka/Ks > 0 implying purifying (negative) selection, neutral evolution and positive selection respectively. Note that this assumes the neutrality of synonymous sites, which is increasingly questioned due codon usage bias and regulatory elements overlapping coding regions (and also ignores mutational biases and biased gene conversion, although these can be accounted for in an evolutionary model, and also assumes no sequencing errors and perfect alignments, which can impact a lot on the results too).


Yes. I haven't found it documented anywhere, but I ran a simple test:

$ bcftools view -Oz --no-version foo.vcf > bar.vcf.gz
$ md5sum bar.vcf.gz
18a2f494710cf170ba79892936015ba8  bar.vcf.gz

$ gunzip bar.vcf.gz
$ bgzip bar.vcf
$ md5sum bar.vcf.gz
18a2f494710cf170ba79892936015ba8  bar.vcf.gz

$ gunzip bar.vcf.gz
$ gzip bar.vcf
$ md5sum bar.vcf.gz
9f5b1cf654966912b78a9cb7175ffc2c  bar.vcf.gz

I passed a vcf file through bcftools view -Oz, checked the md5sum of the resulting file, then uncompressed it and recompressed with bgzip and confirmed that the resulting file had the same md5sum. Finally, I also compressed using normal gzip to be sure that the md5sum would be different, as it was.

So yes, it does indeed look like bcftools creates files that are the same as though you'd passed them through bgzip.


p.adjust can calculate adjusted p-values with different methods. I think you can choose Bonferroni correction and BH correction, for instance. BH correction gives you False Discovery Rate (FDR), which is one of the most commonly used I guess.

Q-value is different and cannot be calculated by p.adjust. The idea is that when you have many p-values to adjust, they form a distribution. If all p-values come from the null hypothesis, the p-value distribution is going to be uniform on [0, 1]. But usually, our p-values are a mixture of the null distribution and alternative distribution. The parameter $\pi_0$ describes the proportion of p-values (or tests) that comes from the null hypothesis. To calculated q-values, the algorithm first estimates $\pi_0$, then calculates the q-value.

Note that if you set $\pi_0=1$, which means that you believe all the p-values come from the null distribution, then q-value is equal to FDR.

I think most people do not know about the meaning and details of q-values and it might be sufficient to use BH procedure in most cases.


Extracting the CIGAR string from a BAM via Python?

If you really do just want the cigar string then it's read.cigarstring. However, I'm not sure what you're trying to gain with the cigar package from Brent. Unless you want to get the string with the masking changed then the tuple you get from pysam is the same as what you get from cigar (with the exception of the numeric operations instead of character operations).


The smooth parameter is an optional keyword and command that changes the smoothing potential. If not provided, the default is 0.5 angstroms by default. Try passing 0 with this keyword in the AutoGrid Parameter File to turn off smoothing potential.

See the AutoDock 4.2.5 documentation for more details.


Your interpretation of pct.1 and pct.2 is misguided. If pct.1 == 1, the given gene was detected in all of the cells. This carries no information about high expression or low expression in the cell, just that it was detected. avg_logFC is the indication of expression levels. So your expectation, that

  "we would instead expect the genes to be expressed at high levels in the cells within the cluster but at low levels in cells not within the cluster"

holds true, as indicated by the positive avg_logFC and significant p-values of the listed RP genes.

Seurat use nautral log, so the FC of RPS6 in cluster 0 vs. all other clusters indicated is 2.718281828459^.55947=1.750.

Seurat FindMarkers() documentation.


Sort reads in BAM file based on presence of specific deletion?

The CIGAR string should tell you if a read has a deletion in it.  So filter the bam files based on that.


What is indel calling and what is its purpose?

There are two types of INDELs: short indels and long indels. Some put the threshold at 50bp; others choose 1000bp. Short and long indels are called differently.

<50bp short indels are called from read alignment directly. Modern indel callers essentially build a multi-alignment of reads and call an indel with enough read supports. Short indels may break open reading frames.

For Illumina reads, long indels are usually called with split alignment or read pair alignment because mainstream aligners are unable to directly align through a long indel. Long indels are a type of structural variation. They may affect gene structures in a more dramatic way (e.g. delete a whole exon; screw up transcription by transposon insertion; create pseudogenes).

Ultimately, both short and long indels are types of genetic variants. Calling them helps to understand how genetics shapes phenotypes.


How to create Seurat object while RNA expression and ADT combined into one matrix

Given the fact that you have RNA and ADT data (probably on the same set of cells), some colnames (which represent individual cells) would be duplicated. Normally duplicated colnames are tolerated with matrices and sparse matrices, however, Seurat apparently does not do so, below is from the ?CreateAssayObject:

Details
Non-unique cell or feature names are not allowed. Please make unique before calling this function.

Possible solution to your problem: Seuart has a dedicated vignette for working with multimodal data and as you would see you will need to initiate your Seurat object with one matrix per assay: RNA and ADT. All you need to do is split your matrix into RNA and ADT, create your Seurat object with RNA data and then add the ADT data with:

seurat_obj_with_rna_only[["ADT"]] <- CreateAssayObject(counts = your_adt_matrix)

For efficiency, Seurat uses sparse matrices so don't forget to convert your data matrices to sparse.


How identifiable are human omics data and how to mitigate their identifying features?

This is one of the major problems with genomic information in todays research. This was highlighted some years ago with the police using publicly available genomic data bases to idenitfy unkown murder suspects.
The full extent of this issue was exemplified in a
paper by Yaniv Erlich from 2018 that should be a good starting point for you.
They claim that giving the current (2018) amount genetic information more than half the searches for a person of european descent will result in a thrid cousin or better match. And these odds will only inprove with time.
Therefore, the genetic information not only can give you a hint at inherited variables like the approximate height of a person but can even yield their name and address.
About your second question, I would think that the amount of host (i.e human) material collected alongside the microbiome sample is rather small. This probably is not sufficient to obtain enough coverage for a profile. But I could be mistaken.


How to represent random sequence elements in SBOL?

This sounds like a job for the CombinatorialDerivation class, which is designed for representing libraries and the relationship between derived designs and the originals.
About CombinatorialDerivation
In SBOL 3, a CombinatorialDerivation takes a Component in its template property to represent the system that you intend to vary:

Some of the Features of the template Component are designated as the variables of the library via VariableFeature child objects of the CombinatorialDerivation, which point at the variables with their variable property.
The values that can be filled in for each variable are then designated via the variant* properties of the VariableFeature.
Finally, the distribution of values is determined by the strategy property of the CombinatorialDerivation and the cardinality properties of each variable.

Representing Random Mutagenesis
In the case of random mutagenesis, if I'm understanding your design correctly, you would represent it with a CombinatorialDerivation as follows:

The template is Component of type DNA with the original sequence. Each site designated for mutagenesis can be designated with a SequenceFeature.
For each of the mutagenesis sites, create a VariableFeature in the CombinatorialDerivation with the corresponding SequenceFeature as its variable and a set of variant values for each of the mutant options set by your primers. The cardinality is probably one (if every site has to be hit by something in your protocol) or zeroOrOne (if sites might or might not be hit).
The strategy for the CombinatorialDerivation is sample, since only a few of the many possibilities will actually be realized.

Representing Samples and Characterized Constructs
Now, when you take it to the lab, you can represent your aliquot with an Implementation whose wasDerivedFrom points to the CombinatorialDerivation, representing the fact that you've got a randomly generated library whose contents are not precisely known.
Finally, once you've got your isolated products represented as Components, there is a recommended best practice to use wasDerivedFrom links to connect the derived Component and its varied Features back to the CombinatorialDerivation and its template. If you make a Collection of your results, the Collection should also link to the CombinatorialDerivation.
This should then capture the relationship between all of the products, their associated data, and the original mutagenesis design.


LD analysis in PLINK based on reference and a SNP list

It looks like the first command is calculating pairwise linkage disequilibrium r2 between your SNPs.

The second command is clumping your dataset based on a p-value over a file of association test results, which I guess is a pre-cursor to computing polygenic risk scores.

Clumping is the process of thinning your SNP list to only keep uncorrelated SNPs so that each representative region has only one marker. In this case you are asking PLINK to keep the marker with the highest p-value.

I'm not sure what your definition of proxy SNPs is but if you are looking for tag SNPs to a SNP of your interest, you should go with command 1 and choose a proper LD r2 threshold.


The answer to this question depends on the aligner used to generate SAM/BAM file. For BWA method that is part of the pipeline mentioned in comments, split reads will be flagged as supplementary alignment (2048). You can modify this behavior by using option -M as described in the BWA-MEM manual page.

So, in case of split alignment, you will end up with multiple aligned segments that have identical query name and they are flagged as 2048. These segments may only have part of the original read sequence, but the complete read sequence is present in at least one of the segments.

Also, there is no unique sign in CIGAR of an individual segment to identify split reads. Each aligned segment of the original read may have long deletion in its CIGAR or aligned pairs. However, without looking at all segments with the read ID it is not clear whether the segment has a novel sequence in it that doesn't align with reference or that part matches different chromosome/coordinates.


How to get fasta alignment file from SAM/BAM file?

I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align inserted sequences. Even if you don't care about inserted sequences, a MSA in fasta is far to big to be practical. Alternatively, by "fasta alignment file", you could mean pairwise alignment, but it is still impractical to output every pairwise alignment in a separate fasta file.

There are tools to convert SAM to BLAST-like format if that is what you want. The following shows one of them. I believe there are other tools easier to use, but I forget what are they.

Anyway, to install:

wget https://raw.githubusercontent.com/lh3/minimap2/master/misc/paftools.js
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` k8

Then you can convert your ALN.sam with:

./k8 paftools.js sam2paf -L ALN.sam | ./k8 paftools.js view -l0 - | less -S

Or generate a SAM stream from ALN.bam and pipe it with:

samtools view -h ALN.bam | ./k8 paftools.js sam2paf -L - | ./k8 paftools.js view -l0 -

Here is an example output:

>MT_orang   16499   0   16025   +   MT_human    16569   576 16569   13773   16095   60  tp:A:P  mm:i:2150   gn:i:172    go:i:101
Ref+:        577 GTTTATGTAGCTTACCTCCT---CAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCA
                 ||||||||||||||  | ||   ||||||||| |||||||||||| | ||||||| |||| | |||||||||||||||||||||||||||||||||||||||||||||||| || |||||||||||||||||||||  |||||||| || ||||| || ||||   ||  || |||| |||||||||
Qry+:          1 GTTTATGTAGCTTA--TTCTATCCAAAGCAATGCACTGAAAATGTCTCGACGGGCCCACA-CGCCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTGAGGTTACACATGCAAGCATCCCCGCCCCAGTGAG-TCGCCCTCCAAGTCACTCTGACTAAGAGGAGCAAGCATCA
//

You can change -l0 to -l80 to print alignments in multiple lines.


You said that you want to compare with other microbial genome data. Have you already annotated your genome? If not, you might consider a pipeline I've developed, I can gladly help if you face some problems using it.

It is called bacannot, it does a complete annotation using Prokka, Kofam, CARD, VFDB, and other resources. It might help you to have a first scratch of your data since it also collapses annotation information in a genome browser


How to represent pathways as mathematical graphs?

Pathways are generally directed graphs, which can be implemented in python using the networkx package. Note that there will be cycles in graphs.

Please note that there are large numbers of pathways analysis packages that already exist (both python and R packages but also just webservers), so I strongly recommend searching to see if something already exists before wasting time and implementing it yourself.


How to you find read statistics for bam files with terminal mismatch?

Could you provide some examples of lines that should be counted? I'm not sure how mismatches are encoded in the version of Bowtie you're using. As far as I remember, Bowtie doesn't distinguish between = and X in CIGARs (and uses M for both), therefore, mismatches must be detected from the MD tag.

You can use SAMsift for filtering the alignments. The following filtering criterion should keep only reads with  read length from 18 to 34 that have an MD tag with 0 as the last character, which must be preceded by a non-number (this should signalise a mismatch at the last position).

samsift \
  -i test1.bam \
  -f '18 <= len(SEQ) <= 34 and MD[-1]=="0" and not MD[-2].isdigit()' \
  -m nonstop-remove

The obtained alignments can be counted by adding | samtools view | wc -l.

Explanation:

-f ... – filtering criterion
18 <= len(SEQ) <= 34 – the allowed read lengths
MD[-1]=="0" – the last character of the MD tag is 0
not MD[-2].isdigit() – the penultimate character of the MD tag is not a digit
-m nonstop-remove – remove alignments without any MD tag (probably unaligned, they would cause an error).


Convert FASTA to FASTQ with dummy quality scores

A simple Biopython solution:
from Bio import SeqIO

for r in SeqIO.parse("myfile.fa", "fasta"):
    r.letter_annotations["solexa_quality"] = [40] * len(r)
    print(r.format("fastq"), end='')

Example:
$ cat << EOF > myfile.fa
>1
ATG
>2
ATGTAGA
EOF
$ python3 myscript.py
@1
ATG
+
III
@2
ATGTAGA
+
IIIIIII


Determining position of side chain hydrogen in glycine residues from coordinates of backbone atoms

Glycine has a single hydrogen atom as its side chain:

All the six bond angles with the CA atom in the middle are about 109°
(C-CA-N, C-CA-HA3, C-CA-HA2, N-CA-HA3, N-CA-HA2 and HA3-CA-HA2 using the CCD naming convention). This defines the rough direction of the hydrogen atoms.

But which of the two hydrogens is the side chain?

If it is a left-handed protein, I'd assume that it's the hydrogen pointing toward the screen, similarly to the side chain in L-alanine:

The images are taken from WikiMedia Commons:
glycine
and L-alanine.


What does the yellow mean in this image from Virus Pathogen Resource?

Sulfur atoms are shown in yellow.

The molecular viewer that you use is JSMol (JMol ported to the web).

Atoms are colored by element: grey C, blue N, red O and yellow S.
If you wonder how other atoms would be colored, see JMol's Default element colors, by periodic table.


Should one be especially conscious of removing low counts, when applying batch correction?

I don't think this has anything to do with batch correction.
The averages of that gene are distinctly different between groups, so it's correct that the gene is statistically significant.  But with one sample driving all the difference, it's likely to be an artifact, or otherwise not fruitful to follow up on.
As geek_y suggested, a filter like "omit genes unless 2 samples have > 10 counts" would drop that gene.


What is the fastest way to get the reverse complement of a DNA sequence in python?

I don't know if it's the fastest, but the following provides an approximately 10x speed up over your functions:

import string
tab = string.maketrans("ACTG", "TGAC")

def reverse_complement_table(seq):
    return seq.translate(tab)[::-1]

The thing with hashing is that it adds a good bit of overhead for a replacement set this small.

For what it's worth, I added that to your code as "with a translation table" and here is what I got on my workstation:

global dict implementation
       1.37599s total,
       363374.8 strings per second
       3.3% increase over baseline
naive implementation
       1.44126s total,
       346919.4 strings per second
       -1.3% increase over baseline
with a translation table
       0.16780s total,
       2979755.6 strings per second
       88.2% increase over baseline

If you need python 3 rather than python 2, then substitute tab = str.maketrans("ACTG", "TGAC") for tab = string.maketrans("ACTG", "TGAC"), since maketrans is now a static method on the str type.

For those wondering, using biopython is slower for this (~50% slower than the naive implementation), presumably due to the overhead of converting the strings to Seq objects. If one were already reading sequences in using biopython, though, I wouldn't be surprised if the performance was much different.


Definition of "seed" in sequence alignment

The seed is the subset of a read used in the first step of an alignment. Many aligners work by a seed-and-extend model, wherein they first find all regions matching the "seed" and then extend the alignment around that allowing mistmatches and indels until it either gives up (and therefore uses a different seed) or finds a sufficiently good alignment.


Mining documents for nucleotide and peptide sequences

I am not sure there is a state of the art for that.
Some possible strategies I can think of:

Applying a dictionary to identify language words and discard these.
Using some context information like identifying some words in the same sentence than the sequence.
I guess the format in papers for sequences will not be the same as for text. Maybe you can find a pattern.

I would use an incremental strategy. Use the simplest algorithm then see what it recovers. If there is garbage, find a strategy to remove some garbage with a common pattern. Repeat until you are happy with the result.
Disclaimer: I have some experience in Natural Language Processing.


They are not referring to density of reads. It is the density of the fold change, the y variable in the plot to left of Fig2b. It tells you the distribution of the foldchange and 92% of the fold change is <0.
See this page for more details on density and histograms
Suppose your data is something like this (purely making this up):
set.seed(555)
fc = rnorm(1000,-0.5,2)
averATAC = -0.5*(fc - 0.5)^2 + 0.2*fc + 3  + rnorm(100,25,10)

You fit a density to it and you get the plot:
par(mfrow=c(1,2))
plot(averATAC,fc)
plot(dens$y,dens$x,xlab = "density",ylab="fc")


Removing data from OrthoFinder output prefix

The latest version solves the issue.
Source : https://github.com/davidemms/OrthoFinder/issues/403
Latest Version: https://github.com/davidemms/OrthoFinder/releases/tag/2.4.0


How does sequence alignment help with protein folding prediction?

Protein sequences from the same family form a sequence alignment. Some positions in the alignment will be conserved, i.e. the same in all sequences. These are often associated with function, e.g. a catalytic residue at position 50.
However there is a more subtle signal in sequence alignments called covariation. Two positions, say residue 30 and residue 60, may appear to change together across the alignment. We infer that they are close in 3D space in the structure:

Initially statistical techniques (DCA, PSICOV etc.) were used to separate direct from indirect (chained) couplings and extract predicted contacts from sequence alignments. Now deep neural networks have improved performance and even allowed the prediction of residue-residue distances and orientations.
These constraints on the structure allow model building, for example using the CNS software. A good set of constraints effectively guarantees the correct fold. The depth and quality of the sequence alignment is the major factor in the accuracy of the predicted constraints.
This approach to protein structure prediction is the current state of the art. Use of covariation from sequence alignments is the only approach that has reliably given accurate models for larger proteins without templates. See this paper and this paper for the best current methods, and tune into the CASP conference this December to see if anyone can improve on them.


Error in colSums(assay(sce)) : 'x' must be an array of at least two dimensions

Ultimately I just implemented it manually with the precursor to SingleCellExperiment matrix (UMI_count which is dgcMatrix type) object in the following way:

# Removing genes with 0 - counts in all of the columns
num_removed = 0
for(i in 1:nrow(UMI_count)) {
  idx = i - num_removed
  if(sum(UMI_count[idx,]) < 1) {
    UMI_count <- UMI_count[-idx,]
    num_removed = num_removed + 1
    print(i)
  }
}

sce <- SingleCellExperiment(list(counts=UMI_count))

I am getting sums of the columns for each row and then filtering out rows with the sum < 1. So, sum function for the slices - rows (numeric type) is working.


I am trying to run nab to create the pdb file but it isn`t working (AMBER)

I have no idea about nab ... however
“$AMBERHOME/bin/nab lr14.nab./a.out”

Looking at your error lr14.nab./a.out:
It suggests you didn't put a space character between lr16.nab and a.out. If you run
$AMBERHOME/bin/nab lr14.nab ./a.out

I suspect it will work because at present it is looking for a file named lr14.nab./a.out which could be interpreted as a directory of lr14.nab. and a input file of a.out or the error suggests it is reading / as part of the file name. a.out of course doesn't exist its your specified output file name

output of ls -l
-rw-r--r-- 1 phys 89 Oct 6 12:58 lr14.nab – Peter 5 hours ago


I'd recommend having a go with either OneCodex (commercial / Freemium) or Centrifuge.

Chris Mason's lab has recently looked at metagenomic classification and OneCodex gave the most accurate result.

Centrifuge uses Burrows-Wheeler lookups, but hasn't yet been compared by the Mason Lab.


This is what I see in my samtools phase (v1.10)

Usage:   samtools phase [options] <in.bam>

samtools phase only takes a bam file as an input. But you are specifying a .fasta file as an input. If you want to specify the reference genome, you should be using the argument "--reference" as described here

      --reference FILE
           Reference sequence FASTA FILE [null]


Which measure should be used in a PCA or RNA-seq data? TPM or counts?

tldr: log transform counts and TPMs, but rlog/vst are preferred

TPM should be log transformed to get more useful results. If you're using DESeq2 already (given the reference to rlog and vst, this seems likely), then please go ahead and use rlog or vst. That will give you more reasonable results than raw counts. If you're stuck with counts for some reason, then first use normalized counts so they're at least a bit more comparable and then log transform them so your highly expressed genes aren't driving everything.

Edit: As an aside, if you know what the batch effect is (e.g., library prep. date), it's sometimes convenient to include that in your model. You can then assess the genes that are actually changed due to that, which is sometimes handy to know (e.g., which genes might be more/less prone to degradation).


From a set operations viewpoint, consider your hg37 regions as a "reference" set and hg38 as a "map" set.

If you have these sets in BED format, BEDOPS lets you do set operations on them with various overlap criteria.

Specifically, you can use BEDOPS bedmap to map hg38 regions to hg37 regions. The bedmap tool lets you assign overlap criteria from as little as one base to as much as full, mutual overlap.

Once you do mapping, you can do operations, like counting the number of regions that meet the overlap criteria between reference and map sets, and calculate percentages.

For example, use --exact to get exactly-matching regions between hg37 and hg38:

$ bedmap --count --exact hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > exact_counts.txt

You can do fractional ("fuzzy") overlaps, where (for example) at least 50% of an hg38 region has to overlap an hg37 region:

$ bedmap --count --fraction-map 0.5 hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > fuzzy_0p5_counts.txt

Getting a percentage out of a "fuzzy" result can be complicated by overlaps between hg38 regions that overlap a reference hg37 region. A more complicated procedure can be implemented with the --echo-map option and some command-line work to implement logic to deal with that case more stringently.

To get a straight-up count of no overlaps, build the inverse — or mask — of your hg38 regions, using Kent tools fetchChromSizes and BEDOPS bedops --difference:

$ fetchChromSizes hg38 | awk '{print $1"\t0\t"$2}' | sort-bed - > hg38.bounds.bed
$ bedops --difference hg38.bounds.bed hg38.bed > hg38.masked.bed

The regions in hg38.masked.bed make up all of the genomic space in hg38, except for your lift-over regions.

Apply --count --fraction-ref 1.0 to count the number of reference (hg37) elements which completely overlap these masked hg38 regions:

$ bedmap --count --fraction-ref 1.0 hg37.bed hg38.masked.bed | awk '{s+=$1}END{print s;}' > no_counts.txt

The documentation for bedmap describes these overlap criteria and operations in more detail.


How to output all sequences with bwa mem, not `*`?

I don't believe that is an option.  Internally we created a patch and submitted it to BWA however that was 4 years ago and it has not been accepted.


The VCF spec clearly states that these are indeed log likelihood values:

  • GL : genotype likelihoods comprised of comma separated floating point log10-scaled likelihoods for all possible
  genotypes given the set of alleles defined in the REF and ALT fields. In presence of the GT field the same
  ploidy is expected and the canonical order is used; without GT field, diploidy is assumed. If A is the allele in
  REF and B,C,... are the alleles as ordered in ALT, the ordering of genotypes for the likelihoods is given by:
  F(j/k) = (k*(k+1)/2)+j. In other words, for biallelic sites the ordering is: AA,AB,BB; for triallelic sites the
  ordering is: AA,AB,BB,AC,BC,CC, etc. For example: GT:GL 0/1:-323.03,-99.29,-802.53 (Floats)

With log-likelihood values, the higher the number, the better the value. So -2 is better than -3.


Download proteomes from NCBI based only on binomial names

The general procedure is explained here:

First step is to search the Taxonomy database with the organism name. Accepted common names usually work at all taxonomic levels. Use the scientific name or formal name if no results are obtained with the common name.

In code:

from Bio import Entrez

Entrez.email = 'user@example.org' # Put your email here

def get_tax_id(species):
    species = species.replace('_', '+').strip()
    search = Entrez.esearch(term=species, db='taxonomy', retmode='xml')
    record = Entrez.read(search)
    return record['IdList'][0]

if __name__ == '__main__':
    organisms = ['lophodermium_seditiosum', 'trichoderma_harzianum']
    taxids = [get_tax_id(organism) for organism in organisms]
    print(taxids)

This gives us the taxids for our species, respectively 128023and 5544. Now we can use thes to get the proteomes. For Lophodermium seditiosum we find two results, but for Trichoderma harzianum there are 65900 results, so you might want to tweak the query to filter the results. For example, if you only want the sequences appearing in PDB you add AND pdb[filter] so that the query becomes txid5544[Organism:exp] AND pdb[filter].

def get_proteomes(taxid):
    query = 'txid{}[Organism:exp]'.format(taxid)
    handle = Entrez.esearch(db='protein', term=query)
    result = Entrez.read(handle)
    if int(result['Count']) <= 2000:
        ids = ','.join(result['IdList'])
        proteomes = Entrez.efetch(db='protein', id=ids, rettype='fasta').read()
        return proteomes
    else:
        raise RuntimeError('Too many results!')

if __name__ == '__main__':
    print(get_proteomes('128023'))

Returns:

>ADI44294.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS

>ADI44293.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS


How to run same command on multiple files?

Another solution for repeating one command with different parameters is to use parallel.

Create a tab-delimited file samples.txt, which contains the sample names and the information for your read groups, e.g.:

Sample1 Lib1
Sample2 Lib2
Sample3 Lib3

Then you can use parallel like this:

parallel --colsep "\t" 'bwa mem -M -R "@RG\tID:{1}\tLB:{2}\tPL:ILLUMINA\tPM:HiSeq2000\tSM:{1}" ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta {1}_* > {1}.sam' :::: samples.txt

For each line in samples.txt a new process is started, replacing {1} and {2} with the values in column 1 and 2 of this line.

fin swimmer


How do I create a filtered gene list using expression medians

This below should work. It returns a list,

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

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

Hope I got it correct!


Are there any Genome in a Bottle-like resources for non-humans (especially for invertebrates)?

In addition I would like to add "Vectorbase", the curator "Dan" is a good guy (I forget his surname right now) and it's a good team. However, please read my rather somber assessment of arthropod genomics first.

We appear to be talking about lower metazoans, e.g. not protists, or yeasts (the latter I suspect are really good for your question). I don't know anything about worms.

Its important to understand that arthropods are several (soccer) divisions lower than mammals in terms of the question you are asking, except for two model organisms.

Flybase, Drosophila melanogaster is the most important species for your purposes. The second most important species within the Class Insecta will be Anopheles gambiae (malaria mosquito).

Personally I would stick to the main players because overall whilst there are a lot of insect genomes (because there are a lot of insects), their genomes remain rough because of the difficulties in assembly and requirement to physical mapping their genomes in comparison to their genome size. Also except for D. melenogaster the (incredible) tardiness in producing a draft genome. For example the idea of PacBio is in its infancy for many entomologists: and this would help them a lot in assembling their genomes correctly. I agree things have picked a lot recently, but I remain a skeptic.

Within the arthropods I would be very surprised if the crustacea, myriapods or Chelicerata (e.g. spiders) were in better condition than hexapods.

If you are interested in protists certainly Plasmodium falciparum (African malaria) has the sort of information you seek (Plasmodb).

Summary
I don't want to discourage you, insects are great, but you need to be selective if you are from a human/mammalian genetics background.


How to replace sequence identifiers in a fasta with OTU IDs from another file?

You can do it with bash quiet easy. First create a backup of your fasta file with
cp otu.fasta otu_backup.fasta

Then you can run:
while IFS=" " read -r otu fastaId rest;
do
sed -E -i "s/>$fastaId .+/>$otu/" otu.fasta; done < ids.txt;
done < OTU_ID.txt

What does the code?
With IFS=" " you declare a whitespace as field seperator and by reading each line of your OTU_ID.txt, you have the two variables otu and fastaId which hold the otu name and the id of the Fasta Header.
Within the while loop, sed is looking each loop for the corresponding id within a RegEx and substitute it with the corresponding otu name. By using the -i option, you overwrite the existing file. This is why you need to create a backup.


How can I output an identity matrix in progressiveMauve?

mauve documentation have a similar command but the usage is --id-matrix=<file>. So very similar to the progressivemauve command (just doesn't include --input). Try it and see if it works.

I think it's a bug for sure as several other people on biostar have complained.


Q: What analyses can I perform with a completely phased genome assembly?

Some population genetic inference methods require the data to be phased because they look for long stretches of matching alleles between individuals, rather than just matching genotypes. Here, 'genotype' data is un-phased, and 'allele' data is phased.

This is because there is a degree of ambiguity when it comes to matching heterozygous positions - if the data is un-phased, you don't know if 2 individuals share mutations on the same chromosome (here, I mean different copies of the same chromosome - like either copy of chromosome 1, for example). This is important, because the fact that the mutation is shared on the same chromosome means that they are much more likely to share recent common ancestry with one another.

Let's say we have genotype data for two individuals (where the number corresponds to the number of non-reference alleles present):

Ind 1

0 ----- 1 ----- 2 ----- 0

Ind 2

0 ----- 1 ----- 2 ----- 0

If were looking at unphased genotype data only, we know the alleles match exactly in the 1st, 3rd and 4th positions, because they are homozgyous (either 0 or 2 copies of the non-referencee allele), but because the second position is heterozygous, we can't be sure which chromosome the non-reference allele is on. This is because the haplotype configurations for individual 1 could be either

0 ----- 0 ----- 1 ----- 0 Chromosome1a
0 ----- 1 ----- 1 ----- 0 Chromosome1b

or

0 ----- 1 ----- 1 ----- 0 Chromosome1a
0 ----- 0 ----- 1 ----- 0 Chromosome1b

i.e., we don't know which chromosome the non-reference allele at position 2 is on.

If we then phased the data, and found that the non-reference allele was on chromosome 1a for Individual 1, and chromosome 1b for Individual 2 there is less evidence of recent shared ancestry than otherwise. I think this is more or less the approach that 23andme takes when they are inferring your relationship to other individuals in the dataset. After writing this, I realised there is a much better description of what I just said on the 23andme website!

For example, the above applies to the ChromoPainter/fineSTRUCTURE algorithms. This means that they are able to detect more subtle genetic differences between individuals than if you were to use non-phased data.

This is a slightly niche example, and there occasions where you need phased data, but this is the one I am most familiar with :). Off the top of my head, methods such as Multiply Sequential Markovian Coalescent (MSMC), RELATE, and Extended Haplotype Homozygosity all require the data to be phased.


Consensus on how to select SNPs for creating a phylogeny when the full dataset is too large

The subset that gives you the tree closest to the one you obtain with the full dataset does not necessarily is the tree nearest to reality... It could be that you get a more realistic result using only partial data for multiple reasons. Your best chance is in my opinion, as it usually is, making lots of tests and contrasting hypotheses.


Why is it necessary to add hydrogen and delete water before protein-ligand docking?

My answer is mainly based on molecular dynamics,but in this context it shouldn't make much of a difference.

The main reason has to do with the force-field-water model used during the simulation. Examples of water models are SPC, TIP3P or TIP4P; same molecule (water) but different topology.

A pdb usually comes from different sources and it can contain molecules that are not standard and therefore not recognised by the model of your simulation. Usually, once the structure is loaded, your system is re-built so that it conforms to the model topology. There are exception though, for example if your water is part of a functional active site. In that case you can keep it, but its structures (topology of the water molecule) has to be recognised by the model.

Hydrogens are usually added to "saturate" and stabilise the bonds.


Variant calling without matched normal sample

If you just want to filter out calls present in dbSNP then use:

java -jar GenomeAnalysisTK.jar \
   -T SelectVariants \
   -R reference.fasta \
   -V patient.vcf \
   --discordance dbSNP.vcf \
   -o patient.filtered.vcf

--discordance will produce calls not present in dbSNP.


Error in cor(exprs(gse), use = "c") : no complete element pairs

The GSE133399 dataset contains both ATAC-seq and RNA-seq data, none of which can be queried with GEOquery which works only for microarray data.
exprs(gse) ## print the expression data returns an empty object as you might have noticed, therefore the errors in the lines that follow this command.
On the bottom of https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133399 at the Supplementary file section you can download FPKM-transformed count data and use these. I would be careful though as FPKM is not ideal for inter-sample comparison. I personally always (if feasable) download the raw data and then process myself, but for a guess on what is going on the FPKMs probably will do.


Should you filter GWAS hits with high standard error?

I am curious about the relationship between your p-value (or effect sizes) and standard error. I would expect the significant signals to have smaller stand error compared to the non-significant, background signals. If this is the case, there is no need to do filter.


Is there an easy way to create a summary of a VCF file (v4.1) with structural variations?

According to the bcftools man page, it is able to produce statistics using the command bcftools stats. Running this myself, the statistics look like what you're asking for:

# This file was produced by bcftools stats (1.2-187-g1a55e45+htslib-1.2.1-256-ga356746) and can be plotted using plot-vcfstats.
# The command line was: bcftools stats  OVLNormalised_STARout_KCCG_called.vcf.gz
#
# Definition of sets:
# ID    [2]id   [3]tab-separated file names
ID  0   OVLNormalised_STARout_KCCG_called.vcf.gz
# SN, Summary numbers:
# SN    [2]id   [3]key  [4]value
SN  0   number of samples:  108
SN  0   number of records:  333
SN  0   number of no-ALTs:  0
SN  0   number of SNPs: 313
SN  0   number of MNPs: 0
SN  0   number of indels:   20
SN  0   number of others:   0
SN  0   number of multiallelic sites:   0
SN  0   number of multiallelic SNP sites:   0
# TSTV, transitions/transversions:
# TSTV  [2]id   [3]ts   [4]tv   [5]ts/tv    [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
TSTV    0   302 11  27.45   302 11  27.45
# SiS, Singleton stats:
...
# IDD, InDel distribution:
# IDD   [2]id   [3]length (deletions negative)  [4]count
IDD 0   -9  1
IDD 0   -2  4
IDD 0   -1  6
IDD 0   1   4
IDD 0   2   1
IDD 0   3   3
IDD 0   4   1
# ST, Substitution types:
# ST    [2]id   [3]type [4]count
ST  0   A>C 2
ST  0   A>G 78
ST  0   A>T 2
ST  0   C>A 5
ST  0   C>G 0
ST  0   C>T 66
ST  0   G>A 67
ST  0   G>C 0
ST  0   G>T 1
ST  0   T>A 1
ST  0   T>C 91
ST  0   T>G 0
# DP, Depth distribution
# DP    [2]id   [3]bin  [4]number of genotypes  [5]fraction of genotypes (%)    [6]number of sites  [7]fraction of sites (%)
DP  0   >500    0   0.000000    333 100.000000


Modifying a SAM header after removing all non-primary reads

Your first command only needs a slight modification to add in -h. This will create a SAM file with a header.

samtools view -h -F 256 input.bam > headered.input.primaryOnly.sam

  -h        Include the header in the output.

If your primary end goal is to create a BAM, then as Pierre has pointed out you can create BAM files directly from samtools view with the -b option:

samtools view -b -F 256 input.bam > input.primaryOnly.bam

And if your end-goal is a sorted BAM file, then pipe the BAM output through samtools sort:

samtools view -b -F 256 input.bam | samtools sort > sorted.input.primaryOnly.bam

So that's fixing what I think your problem is, but what about the specific question that you have asked, i.e. how can a header be added to an existing headerless SAM file?

If you have an existing BAM file, you can use Picard (as you have mentioned), or manually piece the files together with cat, which will have a similar function (although Picard probably includes a few checks to make sure that the header is consistent with the input file):

(samtools view -H input.bam; cat input.primaryOnly.sam) > headered.input.primaryOnly.sam

If the only thing that matters is that the BAM file has correct reference information (which is fairly common), then the -T option can be used with samtools view. Again, -h or -b is required to make sure both the header and the usual SAM fields are displayed:

samtools view -T reference.fa -h input.primaryOnly.sam > headered.input.primaryOnly.sam
samtools view -T reference.fa -b input.primaryOnly.sam > headered.input.primaryOnly.bam


Is there an "official" standard usage of Red and Green on differential gene expression heatmaps?

As others already pointed out, there is no official heatmap color scheme. I am not sure there is an official color scheme for any other types of plots either.
However, heatmaps seem to be far more standardized. Traditionally, heatmaps have been red-black-green as you originally stated. Most heatmaps from the microarray era and early RNA-seq are like that. You can have a heatmap for any type of data, but the majority is for gene expression. A lot of tools and tutorials were written back then, so they have that color scheme. From about 2008 to 2012, most people switched to the red-white-blue (or red-yellow-blue) color scheme. Switching green to blue avoids the color-blindness issue. Switching black to white is better for the ink if you print your papers and generally looks better on white background.
I originally found red as high to be confusing because I was thinking in financial terms (if the market is red, it's down). Now I think of it as a water faucet (red is hot, blue is cold).
Edit: In the last year or so, the viridis palette started to become widely used. It's far from the standard, but could become more popular in the future.


How are residuals used as new phenotype?

The residual is the level of error in a regression model, the lower the residual the better the model. Residuals cannot equate to a phenotype, it is actually the opposite if your regression model is trying to investigate a given genotype/phenotype(s) model. The residual in this case is the amount of infomation which the genotype/phentype does not account for, i.e. there are other factors resulting in your observations which you have not captured in your model.

What I assume is that the regression model you are referring  is  investigating whether the residual is lowered by examining different gene combinations. The objective of a regression model is to lower the residual and I suspect that it is comparing the new residual against the old one. Therefore it is not being used as an input into a regression model but an assessment of its robustness.


Can I run STAR without an annotation file?

I don't think you can use the --quantMode GeneCounts option with no annotations. I think the error is trying to look for an exon file generated from the annotations to do the quantitation on. Remove that and I think it should work, as the manual specifically states that annotations are optional but highly recommended.


P-value correction when evaluating correlation between gene and miRNA expression

I made the same question in CrossValidated forum and got an excellent answer!
The important part:

You need to correct for all of the comparisons you are doing. So if
that's 300,000,000 comparisons you need to correct for that many
multiple comparisons.

For more information check the answer in the link above


What unit is the length of a protein measured in?

As mentioned by others, it's the number of amino acids. Alternative measurements could be in Angstroms or similar, but that would require at least an approximate structure.


The best would be Biopython using AlignIO, given they are all in a single file otherise SeqIO. The alignment approach,
from Bio import AlignIO
align = AlignIO.read("/home/path/myseq.fa", "fasta")
print(align)
with open("/home/path/example.nex", "w") as handle:
    SeqIO.write(align, handle, "nexus")

(not debugged)
MEGAX or ClustalX which are nice GUI's that will also work. If you've not got Biopython installed then ClustalX is probably the easiest.


The Bioinformatics Stronghold at Rosalind is an excellent resource for building your understanding of many concepts and algorithms you'll see frequently in bioinformatics and genomics.

I'd also say that, if you can make it work, an internship where you can "get your hands dirty" with real problems is a very valuable experience. This option may not be available to everyone due to a variety of factors, but when it is an option it provides a tremendous opportunity to build experience with little risk to the intern or the employer.


RNAseq: Z score, Intensity, and Resources

It depends on what test or analysis you want to do, whether you need intensities (expression values) or z-scores.

If you want to do statistical analysis, such as finding differentially expressed genes between groups of patients (e.g., with limma), you don't want to use z-scores. But you use normalized intensities (for microarrays) or start with raw counts for RNAseq. Then follow the user guide (limma has a great user guide, also for starters).

For visualization in heatmaps or for other clustering (e.g., k-means, fuzzy) it is useful to use z-scores. Z-scores are a form of transformation (scaling), where every genes is sort of "reset" to the mean of all samples, using also the standard deviation. If you want to know exactly what a z-score is, a simple google search can tell you the details.

In R you can use the scale function for z-score transformation. Be aware that the function works on columns though. Which means you have to transpose your matrix first if genes are in the rows, and then transpose them back after transformation.


How do you find differences when your metagenome sample size is small?

I think there are several ways to figure this out. I wouldn't worry too much about normality as most of the methods are either nonparametric or pretty robust to violations.

Since you are using KEGG, I would suggest working at the pathway or module level of annotation (if your features are e.g. KEGG Orthology groups). In principle, you can then combine information across features as well as across samples. Simplest I can think of, without even looking at the literature, would be to use permutations (shuffle row and column labels) to look for modules different in abundance between groups. (Many existing methods are doing some more sophisticated version of this.)
I think that LEfSe in theory should work quite well here, it's been a while since I thought about it but I believe this is the problem that it is designed to solve, assuming that I understand your question.
There are a lot of methods that are previously described for doing this kind of work. I would suggest reading some reviews of methods, for example the section "Multivariate differential abundance testing" in this article, or here, looking for some inspiration. To some extent it depends on which questions you are interested in answering and what exactly your features represent.


How to predict the distance between two residues in a protein sequence

Developments in this area over the last year are worth pointing out. All use residue-residue covariation information from a protein family so require a multiple sequence alignment to work, though it sounds like you have that.

The most accurate way to do this would be to use pre-trained machine learning methods, which would not require you to train your own version. Examples include:

DMPfold from our lab, which is freely-available open source.
The RaptorX method from the Xu lab.
The state of the art is AlphaFold from DeepMind, though this has not been published yet.

These predict the probabilities of a residue-residue distance being in various distance 'bins', e.g. the probability of being 5 - 6 Angstrom apart. They don't predict beyond ~20 Angstrom.

If you specifically wanted a statistical, non-trained method, then consider:

PSICOV from our lab.
CCMPred.
EVfold.

All of which are available in some form. The problem here is that they all predict binary contacts at 8 Angstrom, whereas it sounds like you need more fine-grained distances.

In a general sense, you can use any 3D structure prediction method and then read the distances back from the predicted structure.


Proper use of BWA MEM on multiplexed GBS sample

The first thing to do is to group the FastQ files together following this post. Once you have your reads grouped, you could trim your files before mapping. Then you could perform your BWA alignment. Since you mentioned MEM, I suppose your read length should be between 70-100bp as suggested by the author:
It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM.
The first algorithm is designed for Illumina sequence reads up to 100bp, while the rest two for
longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar features such
as long-read support and split alignment, but BWA-MEM, which is the latest,
is generally recommended for high-quality queries as it is faster and more accurate.
BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.

Then you could run your BWA MEM, supposing you have 10 CPUs available, like so:
bwa mem \                                  # bwa mem algorithm
-M \                                       # mark shorter split hits as secondary (for Picard compatibility).
-t 10 \                                    # CPUs thread
reference.fasta \                          # reference genome to be aligned to
< sample_r1.fastq \                        # fastQ read 1 (can be .gz)
< sample_r2.fastq | \                      # optional FastQ read 2 (can be .gz)
samtools view -bS -@ 10 > sample.bam       # pipe to samtools to convert SAM to BAM on the fly.

Hope this helps !


How to filter out cross alignments from a BED file?

According to the SAM specification, the 3rd field of a SAM line (RNAME) is:

  RNAME: Reference sequence NAME of the alignment. If @SQ header lines
  are present, RNAME (if not ‘*’) must be present in one of the SQ-SN
  tag. An unmapped segment without coordinate has a ‘*’ at this field.
  However, an unmapped segment may also have an ordinary coordinate such
  that it can be placed at a desired position after sorting. If RNAME is
  ‘*’, no assumptions can be made about POS and CIGAR.

And the 7th field is (emphasis mine, missing "to" theirs):

  RNEXT: Reference sequence name of the primary alignment of the NEXT
  read in the template. For the last read, the next read is the first
  read in the template. If @SQ header lines are present, RNEXT (if not
  ‘*’ or ‘=’) must be present in one of the SQ-SN tag. This field is set
  as ‘*’ when the information is unavailable, and set as ‘=’ if RNEXT is
  identical RNAME. If not ‘=’ and the next read in the template has one
  primary mapping (see also bit 0x100 in FLAG), this field is identical
  to RNAME at the primary line of the next read. If RNEXT is ‘*’, no
  assumptions can be made on PNEXT and bit 0x20

So, you want to remove those lines whose 7th field isn't = and, just in case, those lines whose 7th field isn't = and isn't the same as the 3rd field. You can therefore use something like this:

samtools view -L test.bed test.bam | awk '$7=="=" || $3==$7

And, to save as a bam file again:

samtools view -L test.bed test.bam | awk '$7=="=" && $3==$7 |
    samtolls view -b > fixed.bam

On a separate note, there's very rarely a need to chain multiple grep commands like that. You can just use \| (or | with the -E or -P options) to separate them. Something like:

samtools view -L test.bed test.bam | grep -v 'chrx\|chr2\|chr10\|chrN'

Or

samtools view -L test.bed test.bam | grep -Ev 'chrx|chr2|chr10|chrN'


If I read quickly thru the paper, I see you refer to Fig 5D. In the methods they describe that they use Fluidigm (i.e. RT-qPCR) for gene expression, not RNA-seq. The data link you give is for small-RNA-seq.


What's the best way to download data from the SRA? Is it really this slow?

Proximity to NCBI may not necessarily give you the fastest transfer speed. AWS may be deliberately throttling the Internet connection to limit the likelihood that people will use it for undesirable things. There's a chance that a home network might be faster, but you're likely to get the fastest connection to NCBI by using an academic system that is linked to NCBI via a research network.

Another possibility is using Aspera for downloads. This is unlikely to help if bandwidth is being throttled, but it might help if there's a bit of congestion through the regular methods:

https://www.ncbi.nlm.nih.gov/public/

NCBI also has an online book about best practises for downloading data from their servers.

On a related note, just in case someone sees this and EBI/ENA is an option, there's a great guide for how to do file transfer using Aspera on the EBI web site:

https://www.ebi.ac.uk/ena/browse/read-download#downloading_files_aspera

  Your command should look similar to this on Unix:

ascp -QT -l 300m -i <aspera connect installation directory>/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:<file or files to download> <download location>

In my case, I've just started downloading some files from a MinION sequencing run. The estimated completion time via standard FTP was 12 hours for about 32GB of data; ascp has reduced that estimated download time to about an hour. Here's the command I used for downloading:

ascp -QT -l 300m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/ERA932/ERA932268/oxfordnanopore_native/20160804_Mock.tar.gz .


If I modify a PDB file with a specific mutation, how to minimize energy?

Yes you can modify the reference PDB file and look for the changes and for this purpose you need visualizers. One among these is Chimera. You can easily carry out the energy minimization steps using Chimera, first mutate your PDB file manually and then provide it as an input for Chimera and the steps for doing energy minimization are fairly simple as given in Chimera docs here.


Output when compute GC content is different between using user input and using string

Is it cheating to use Biopython(?), will a solution might be helpful anyway. This assumes you have your sequences stored in the seq.fa FASTA file:

from Bio import SeqIO
from Bio.SeqUtils import GC

print(max(((r.id, GC(r.seq)) for r in SeqIO.parse("seq.fa", "fasta")),
          key=lambda x: x[1]))

Output:

('Rosalind_0808', 60.91954022988506)


Output file names are not known beforehand

Could you give an example of names of input files and output files of normalize.py ?

It seems to me that your case could be similar to the one described about the dynamic output:

https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#dynamic-files

Here you would probably use it for the output of normalization and for the input of all (with a wildcard different from hdf5_file).

Disclaimer: I haven't had the opportunity to use this feature so far.


Based on the comments of the question I thought to put together this answer using GNU coreutils which should exist on pretty much all systems.

The idea is to paste the 2 files together with paste (assuming they are in the same order and each line correspond to the same line in the other file).

Since only the sequence and quality of both files should be kept, the lines starting with @ and + need to be made empty for one of the files. This can be achived with sed, e.g. sed 's/^@.*//' empties all lines starting with @ but keeps the empty line.

To modify the second file on the fly, one can using process substitution <().
Everything strung together should be something like this:

file1=File_A.fastq
file2=File_B.fastq
paste -d "" $file1 <(cat $file2 | sed 's/^@.*//' | sed 's/^+.*//') > target.fastq

I haven't tested it thoroughly, but could be a good starting point.

If you have a lot of files to paste, you might however want to process them first (not on the fly) using:

for file in $(ls *.fastq); do
  cat $file | sed 's/^@.*//' | sed 's/^+.*//' > ${file%.fastq}_cleaned.fastq
done

This creates File_A_mod.fastq from File_A.fastq etc.
Best to move them in a separate directory (make sure the first one is the original instead of the _cleaned one since it serves as the base).
You should have the following in the folder: File_A.fastq File_B_cleaned.fastq File_C_cleaned.fastq ....
Those can be than pasted with paste -d "" *.fastq > target.fastq.


ATAC-seq peak annotation and downstream analysis

The way we do this (and I think its quite common?), is to merge all the files of the same condition (taking the same number of reads from each), and then calling peaks on the merged sample - so you would have two peak sets - one from normal and one from disease. You then do the union of the two peak sets.

This makes annotation much easier.

For downstream quantitative analysis, you can then quantitate each individual sample against the merged peak set.


Is there a publicly available tumor-normal sample?

Please take a look at

An open access pilot freely sharing cancer genomic data from participants in Texas

Although I have not worked on this paper but author claimed that both tumor and normal data is publicly available. In case if you still need something more, than you can work on virtual normal data. And there are several papers which provide machnine learning techniques to overcome the deficiency of normal data (not much efficient but still they give considerable results). Hope this information will help you. Thanks


tabix indexing a GTEx variant file- ':' delimited

Tabix files need a name and a genomic coordinate at least. Tabix can be (ab)used by creating a fake genomic coordinate and then it just indexes names but not sure I would recommend this usage. What is your actual expected application of this?