Understanding SNP coding for association analysis

Humans are diploid. This means they have two copies of each locus (with the exception of sex chromosomes in males). Your two alleles for each individual represent the two copies of this variant locus.


How do I get the adjacency matrix from gene similarity network?

The STRING download page is here. I am not sure that they provide any data in the form of an adjacency matrix, as the number of nodes (rows/columns) would be quite large, and for this reason they are likely to represent the data as lists of links, for compactness.

Their data will likely have to be cleaned for any purpose you are interested in. You will have to make some decisions on how to filter the data. For example, as suggested in the comments, you may want to filter on the quality of the evidence supporting interactions, or based on species of the protein, or other things.

It's possible that you may have better luck going with another database that has done some of these steps already. For example, there are relatively compact human interactome datasets here, or yeast interactome datasets here. In this case you will still need to convert from a list of links to an adjacency matrix, however this is not very difficult if you can use a scripting language.

Notably, you may have to figure out how to distinguish missing data from tested non-interactions in these cases (or decide that missing data is ok). In most cases people will not have exhaustively tested every cell of the all-by-all interactome adjacency matrix.

Hope that helps.


What causes the difference in total length of assembled contigs and scaffolds in SOAPdenovo2?

Something else that could be happening is that contigs that are being collapsed into "heterozygous" groups. This would be a particular problem when a genome has a substantial amount of repeated sequence. Digging deep into the supplementary information of the SOAPdenovo2 paper, I've found the following information:

  In SOAPdenovo2, heterozygous contig pairs are recognized by utilizing the information of contig depth and the locality of contig. The recognized heterozygous contig pairs should obey the following rules: 1) the similarity between contigs should be high enough, for example, ≥ 95%; 2) the depth of both contigs should be near half of the average depth or all contigs, complying Poisson distribution; 3) the two contigs should be located adjacently in a scaffold and have no relationship to each other inferred by paired-end reads information. The normal contigs neighboring the heterozygous regions, if they exist, could be connected to both of the heterozygous contig pairs (H1 and H2). Only the contig with relatively higher depth in a heterozygous contig pair were kept for scaffolding. The method reduces the influence of genome heterozygosity on final scaffold length. All heterozygous contig pairs were outputted to a file to facilitate further analysis. However, the trade-off of this method is that it might incorrectly remove paralogous contigs. This problem could be relieved by a gap-filling procedure while the removed copy of paralogous contigs would be represented by gaps during scaffolding.


Does the kinship and inbreeding coefficients depend on population frequency of an allele?

The definition as set by Kenneth Lange's book is indeed quite vague and thus unrigorous. It is a probability not conditioned on observing any particular allele at a locus thus not on the population frequency of that allele, but conditioned on the following where we focus on one particular locus. It is rigorously defined as follows
A graph of lineage contains nodes, which we call persons, of a particular locus. Each node, or person, contains exactly two genes. We draw one and only one directed edge (parental relationship) from a gene of a parent person to the gene of a child person if the latter is inherited from the former.  This is a directed acyclic graph (DAG). Each gene can be connected by no more than $1$ edge directed to (as opposed to away from) it. The two genes in one person cannot be connected by edges. Two persons are said to be connected if there is an edge pointing from a gene of one node to a gene of the other node. Two persons cannot be connected by more than one edge. For a pair of connected persons and their directed edge between them, the edge is equally probable to connect to one gene as the other gene of a person, i.e. the probability is $\frac12$. A pair of genes (persons) are said to be connected if there exists an undirected (simple) path from one gene (person) to the other. Two thus connected genes are called identical by descent.
In this setting, the DAG with respect to the persons is given and deterministic whilst any directed edge in the DAG assigned to two persons are random with independent probability $\frac12$ with respect to the genes within each person.
For a pair of persons on a lineage DAG that is connected, the probability of one gene equally probably chosen amongst the two genes of one person being connected, or identical by descent, to one gene equally probably chosen amongst the two genes of the other person is the coefficient of kinship. It can be proven that the probability is $\displaystyle\sum_{p\in P} \frac1{2^{n(p)}}$ where $P$ is the set of all distinct (simple) paths connecting one person (not gene) to to the other and $n(p)$ is the number of persons on person connecting path $p$.
As an additional assumption, the genes that are connected have to be the same and thus assume the same allele.

With regard to the conditional probability example in the question, the coefficient of kinship is not adequate to be used to compute the desired probability. Part of the reason is the question inquires about the existence of an allele rather than picking one gene with 1/2 of the probability. We need some more fundamental probability to compute that. We will demonstrate that as follows.

We will compute the probability of the two siblings having both allele $a$ as follows.

For a pair of full siblings, they can have exactly $0, 1, 2$ identical by descent genes. Let the respective probabilities be $r_0, r_1, r_2$. Notice we are not picking genes of a person at random but talking about the existence of particular types of genes in a person. For two full siblings, some thoughts lead to $r_0=\frac14, r_1=\frac12, r_2=\frac14$.
By symmetry and without the loss of generality, we fix the connecting genes in the two siblings of the respective edges. We only need to consider the connections of the edges relative to the genes of the parents.
(1) Conditioned on there is no identical by descent genes, the sibling genes can independently assume any allele. The probability of having the configuration of $ab|ab$ (not considering gene ordering) is $2^2a^2b^2$ considering the distinct ordering of the alleles, that of $aa|ab$ is $\binom41 a^3b$, and that of $aa|aa$ is $a^4$. The total probability is thus
$$Pr_0:=\frac14(4a^2b^2+4a^3b+a^4).$$
(2) When there is exactly $1$ pair of identical by descent genes, the two identical by descent genes can only assume the same allele, while the other genes can assume alleles independently. When the identical by descent genes assume allele $a$ the probability of the two siblings both having allele $a$ is $a$. When the two identical by descent genes assume allele $b$ the desired probability is then $ba^2$. The total probability is then
$$Pr_1:=\frac12(a+ba^2).$$
(3) When there are exactly $2$ pairs of identical by descent genes, each pair can only assume the same allele but independently of each other. The only way that the two siblings not to both have allele $a$ is for the two identical by descent pairs to assume allele $b$. The total desired probability for the two siblings to both have allele $a$ is then
$$Pr_2:=\frac14(1-b^2)$$

Now we compute the probability of any one of the siblings having allele $a$. Each gene can of course assume alleles independently. So the probability is simply
$$Pr_3:=1-b^2.$$

Therefor the desired probability of both siblings having allele $a$ given that one of them has allele $a$ is
$$Pr(a):=\frac{Pr_0+Pr_1+Pr_2}{Pr_3}=\frac{\frac14a(a^3-6a^2+5a+4)}{a(2-a)}$$
since $b=1-a$. Some tedious algebra which can be carried out by, say, Mathematica shows this is identical to the expression in the question in terms of $r=\frac{a}{1-a}$.
$$Pr'(a)=\frac12\Big(2-\frac1{(2-a)^2}-a\Big).$$
It is easy to see $2Pr'(a)-2$ strictly decreases for $a\in[0,2)$ as the last two terms decreases in that interval. $Pr'(1)=0$ and thus $Pr'(a)>0,\ \forall a\in[0,1)$. $Pr(a)$ thus strictly increases for $a\in[0,1]$ and $\frac12=Pr(0)\le Pr(a)\le Pr(1)=1$.


Plink 1.9 equivalent of --score sum in Plink 2.0

plink2 --score cols=+scoresums

See https://www.cog-genomics.org/plink/2.0/general_usage#colset and https://www.cog-genomics.org/plink/2.0/score for details.


Filter reads belonging to unique sequences with threshold

Your algorithm generally works, you just had one minor issue, and a few stylistic changes that could be made.

I'd adjust your script to the following:

for fq in fastq1/*.fq.gz; do
    basename=$(basename "$fq" .fq.gz)
    zcat "$fq" | awk "/^@/{getline; print}" | sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt
 zcat | grep --no-group-separator -B 1 -A 2 -f tmp.txt "$fq" | pigz -p 10 > "fastq2/$basename.fastq.gz"
done

rm tmp.txt

Instead of

for fq in $(ls fastq1/*.fq.gz)

use this

for fq in fastq1/*.fastq.gz

Iterating over ls is generally not ideal, and can give unexpected results. Better to just rely on bash globbing.

Purely stylistic, but I'd simplify and use basename to remove the suffix, instead of cut.
Rather than removing the separators that grep will introduce by filtering them out, use the following option to not create them:
grep --no-group-separator
And, minor bug, as others have mentioned, you should use zgrep or zcat | grep since your files are compressed.

Oh - and quote your variables.

Finally, while this does work, it may not be the most efficient method. The file is being decompressed twice, and may be better to decompress once to a temporary file that can be used and then deleted.


True bam file encoding and viewing as text

Why not work only with sam files and we need to convert them to bam?

SAM is for human only, a binary format like BAM is smaller and parsed much faster by a program (and we want speed and optimize spaces for those huge files, don't we ?).

Furthermore, bam files are compressed using a  method that allow them to be accessed using random-access, that is how you can query a sub-part of bam file for a given genomic location.

you can always go back to SAM using samtools.


Why Illumina if PacBio provides longer and better reads?

There are so many reasons why one might want to prefer Illumina over PacBio (also note that it's a false dichotomy, at least Oxford Nanopore is a competitive sequencing platform):

The first (IMHO and the most common reason) is still the cost of both sequencing and the instruments. Illumina can sequence a Gbp of data for \$7 - \$93. PacBio sequencing is according the same webpage \$115 per Gbp, however at our sequencing center it's ~$200. Though ONT might have already a cheaper solution. Edit, I just found a google sheat with prices that seems to be frequently updated, seems that the ratio still holds Illumina short reads ~10x cheaper than PacBio.
RNA-seq (i.e. analysis of a gene expression) is not possible with PacBio due to preferential sequencing of smaller fragments; shorter genes would always be shown to be more expressed. To be clear, it's possible to sequence RNA with PacBio (the keyword is iso-seq), but the analysis of gene expression is problematic.
It's way easier to extract fragmented DNA (concerns small non-model organisms; although recently a single mosquito was sequenced, so we can expect a further improvement)
other sequencing techniques as RAD-seq that allow genotyping with very little effort and cost, I have never seen anybody even considering using long reads for such genotyping
Genome profiling (assembly-less genome evaluation) based on kmer spectra analysis is not possible with PacBio data due to higher error rates. Conflict of interest: I am a developer of one of the tools for genome profiling (smudgeplot)
DNA in Formalin-Fixed Paraffin-Embedded samples is fragmented anyway (~300bp), therefore, there is no point in sequencing them with more expensive long-read technology (contributed by @mRoten)
I bet there will be a plenty other applications I am not aware of


What assembler is appropriate for High-Fidelity PacBio reads

As recently tweeted by A. Phillippy, Canu v1.9 now supports HiFi reads.

Therefore you should just read their manual and look for the implementation

:)


As you can see on this blog post of the author, the result is ordered by pvalue, but it can coincidentally sort by setSize. The setSize means that that term has so many genes on under that category.


Convert or extract Seurat object as the input of FateID

Seurat objects have various slots where different types and forms of data are stored. For FateID you would need an expression matrix (with cells as columns and genes as rows) and this can be extracted via:

seurat_object@assays$RNA@data

Here RNA is the assay (you can have multiple assays per sample, RNA, protein, ...) and data is normalized data. You can extract counts or scaled data instead.

If you are using Seurat v2, drop the @assays$RNA part.


How to get all chemical compounds that interact with a given protein target?

The easiest solution to the problem of mapping chemical compounds to protein IDs is to use the ChEMBL compound ID to UniProt mapping.

Doing this with PubChem is not so straight forward in part because its downloadable data format is a hideous obfuscated and badly documented mess split over multiple compressed files which in turn contain multiple compressed files...... For a publicly funded project they should really do better.


Is there a simple command for outputting a tab delimited columns?

I have no idea what do you want to do with fasta files reformatted like this, but I think something like this should do the job (untested):
from Bio import SeqIO
import sys

for olig_fasta in SeqIO.parse(sys.argv[1], "fasta"):

  name = olig_fasta.name
  seq = str(olig_fasta.seq)
  seq_len = str(len(olig_fasta.seq))

  print("\t".join([name, seq_len, seq]))


Determining significance of a variable in a glm model

There is no better method, it's a matter of what you want to test or what is your question.
Using the anova() or aov(), test the terms collectively. For example, in your example with Substrate, the null hypothesis is that the coefficients are all zero, meaning cellulose =0, mannan =0 , ....
If the question is, "do the isolate samples have a higher value than origin samples?", then you can use your first model, where free is set as the reference and you test whether the effect of isolate is non-zero. Likewise you can do this for substrate and set of them as your reference. You can also do other pairwise comparisons using this model.
If the question is, "does origin have a significant effect on value, after controlling for substrate?", then you can use your second model.


how to find the bound form of an enzyme structure?

Tried looking for an explicit database? i.e.

ComSin: database of protein structures in bound (complex) and unbound (single) states in relation to their intrinsic disorder: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2808974/ or
LigASite: a database of biologically relevant binding sites in proteins with known apo-structures https://www.ncbi.nlm.nih.gov/pubmed/17933762

Else the PDB itself has a very powerful advanced search that should allow you to do this: http://www.rcsb.org/pdb/search/advSearch.do


Subsetting from seurat object based on orig.ident?

If you are going to use idents like that, make sure that you have told the software what your default ident category is.
This works for me, with the metadata column being called "group", and "endo" being one possible group there.
Idents(combined.all) <- "group"
endo_subset <- subset(combined.all, idents = c("endo"))

I think this is basically what you did, but I think this looks a little nicer


Retrieve RNA sequencing data for human p53 colon cancer cell lines

I suggest you to try difference resources, including ArrayExpress.
A simple search at ArrayExpress returns multiple experiments which might be of your interest. Just a pointer, ArrayExpress is also synced with GEO.


How could I match atom orders between a .mol2 and a .pdb?

Atom reordering is a common problem in compchemistry.

Rdkit (a python package) can do this, but it is limited by the formats it can read and mol2 files are a bit hit or miss. It works really well with SMILES, SMARTS and mol (sdf) files. But the writing may cause problems with Brookhaven pdb and mol2 files.
So the formats you have are both problematic. One is a pdbqt (a PDB file with partial charges (Q, generally Gasteiger calculated) and atom type (T) and the other a mol2. Let's pretend first that they are not.

Also the hydrogens might be made implicit depending on your mol2 file.

from rdkit import Chem
ref = Chem.MolFromMol2File('ligand.mol2')
target = Chem.MolFromPDBFile('ligand.pdb')
atoms = target.GetSubstructMatches(ref) #: List[int]
new = Chem.RenumberAtoms(target, atoms)
Chem.MolToMolFile(new, 'neoligand.mol')

Extra steps would be required for the protons and if your starting ligand was not from a file, but a SMILES.

mol = Chem.AddHs(mol) #protonate explicitly
Chem.GetSSSR(mol) #resonance fixing
AllChem.EmbedMolecule(mol) #initialise for 3d.
AllChem.UFFOptimizeMolecule(mol, maxIters=2000)

So this is all nice, but the formats may not work. Luckily openbabel can deal with this problem.

obabel -i mol2 input.mol2 -o pdbqt -O output.pdbqt

However, one thing to keep an eye out is what will be recalculated. Say, sdf files do not have partial charges or what PDB calls "atom names" and Rdkit calls "atom labels" (About the latter, I actually wrote a script to relabel atom labels in pdbs from mol files), while mol2 and PDBQT have.


Find rsIDs for GRCh37 SNPs and rsIDs for GRCh38 SNPs and compare for overlap

it's not clear to me that it would support both versions of genome
  builds.

You can specify what version of the database you are querying.

I'd start with biomart on the ensembl website, see if you can get it to do what you want, then work on getting biomaRt to do what you want.


Visualisation of long read RNA-Seq splicing

DNASTAR's software is for purchase, but high quality. GenVision Pro does genomic visualization, including Sashimi plots.

Edit: not sure why this answer is being downvoted, unless it's because the software isn't free. OP has tried IGV and SeqMonk, I mentioned an alternative he might not have heard of.

Here is a video demonstrating the use of Sashimi plots in GenVision Pro:

http://www.dnastar.com/t-support-videos.aspx?video=YJvcERoSIsg


Align the HIV-1 Gag protein with the Gag protein of Visna virus

Both are lentiviruses and env is a fast evolving gene. What you are looking at here is the location of indels and this will vary between alignment algorithms.
There will always be an element of error in attempting to align over this genetic distance. It depends what you end goal is, if it is to maximise homology than okay, however optical enhancement of alignments are common if your downstream analysis is phylogenetics. If it is protein structure homology style modelling, where you have the protein structure of HIV env experimentally determined mulitple times you want to use it to help estimate that of the cat lentivirus then optimising a scoring metric is okay.
What you have to understand is there is no 'right' method, in phylogenetics this doesn't matter as much because indels are excluded in the subsequent calculation, so it is less sensative to indel positioning errors and optical enhancement is used. You can also simply exclude regions that appear 'garbled'. For protein homology modeling this is a big problem, because the indels will make drastic changes to the structure prediction. So the only approach here in my opinion is to produce a separate structural prediction for every variation in indel alignment you have obtained. In the latter context getting the alignment right is everything and I suspect thats the output here, because its a classic surface antigen to understand, which will readily crystalise.


How to find a homolog that has a PDB structure available

You could use HMM-based comparison, for example with HHpred. This is generally the best way to find structural templates because structure is more conserved than sequence, and comparing HMMs will turn up hits of the same structure but minimal sequence overlap.

If HHpred turns up nothing you are in the region of de novo modelling, which is a different problem. However, it has advanced in recent years to the point where it may be useful to you. For example, you could try our method DMPfold.


Biopython Genbank writer not splitting long lines

I read the Biopython Genbank record source code carefully and realized I wasn’t using the Qualifier class correctly.

After reading the source code again, I realized I was using the data structure incorrectly.  I need to use the Qualifier class to add a key and feature.  This is formatted correctly and solves the problem.

feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
feature.key = "misc_feature"
#need to add 1 to start list since csv counting is off
feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i])
liability_qualifier = Qualifier()
liability_qualifier.key = '/note="'
liability_qualifier.value = '\"{}\"'.format(liability).split(":")[1]
feature.qualifiers.append(liability_qualifier)


Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your samtools view -H code. With samtools idxstats file.bam you'll get the reads per chromosome (or per nucleotide sequence in your case).

For use in featureCounts my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).

If your data is properly aligned to hg19 your output should be something like this:

samtools view -H aligned2hg19.bam

@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
@SQ SN:chr8 LN:146364022
@SQ SN:chr9 LN:141213431
@SQ SN:chr10    LN:135534747
@SQ SN:chr11    LN:135006516
@SQ SN:chr12    LN:133851895
@SQ SN:chr13    LN:115169878
@SQ SN:chr14    LN:107349540
@SQ SN:chr15    LN:102531392
@SQ SN:chr16    LN:90354753
@SQ SN:chr17    LN:81195210
@SQ SN:chr18    LN:78077248
@SQ SN:chr19    LN:59128983
@SQ SN:chr20    LN:63025520
@SQ SN:chr21    LN:48129895
@SQ SN:chr22    LN:51304566
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566


I got the answer from Rory Stark. the samples should have a unique ID. (https://support.bioconductor.org/p/132502/#132514)


Too many requests error when querying NCBI

The execution of your EDirect commands may be limited unless you are using an API Key, which you can create from your MyNCBI account. Your API Key can be used in one of several different ways.
In an argument to an EDirect command:
esearch -db nuccore -query 'some query' -api_key 12345

In a UNIX environment:
You can either run this just before you are running some edirect commands or add the following line to your .bashrc file.
export NCBI_API_KEY=12345

In the environment, from a Perl script:
$ENV{NCBI_API_KEY} = "12345"


Visualize similar interactions in two homologous networks

An easy way to visualize this is to make a tab-delimited table that contains edge information for both species and both networks, color the edges depending on the species, and make the thickness and transparency dependent on the interaction strength.

For example, say you have a list of orthologs for the two species and just assign them a number or a name. Each edge has a "source" and a "target". The edge weight is the interaction strength.

(Note do not copy-paste the below table as the mixed tabs/spaces may make the graph topology incorrect). In my toy dataset the data look like this:

source  target  variable    value
gene1   gene2   speciesA    3
gene1   gene3   speciesA    42
gene2   gene3   speciesA    3
gene1   gene2   speciesB    9
gene2   gene3   speciesB    5
gene2   gene4   speciesB    41
gene3   gene4   speciesB    32

Import your network to cytoscape as a network from a file:

Then go to the style tab:

... and make sure edge is selected.

Then you can adjust the transparency, edge color, and edge stroke settings, like so using edge stroke color:

:

... and transparency panels as examples.

Don't forget to save your style when it is complete!

When you have adjusted these settings you will have something that looks like this:


Official documentation for textual BLAST output?

I think gringer is correct to say this is because of the BLOSUM62 matrix and a ‘+’ symbol inidcating a positive score for a mismatch but here is a little more explanation.

BLASTP looks for (nearly) identical short protein matches pairwise, called seeds, and then extends the seeds in each direction one residue at a time using arbitrary gap penalties for indels and a matrix for substitutions. Based on some thresholds and more, the highest scoring alignments are kept and ultimately reported.

The default choice of matrix for BLASTP is BLOSUM62, which looks like this. Obviously matches have high positive scores, but also mismatches have varying positive or negative scores.

If you look up your changes:

H -> Q; ; score 0

E -> K; +; score 1

Note that this matrix is symmetric so the scores for H -> Q == Q -> H and E -> K == K -> E.

In terms of why E -> K has a positive score, well not all amino acid changes are viewed equally. BLOSUM matrices are not based on any biological theory but instead use the log-odds ratio, a likelihood approach, that estimates from real sequence alignments how frequently different amino acid changes occur. The more often E -> K substitutions are observed in these alignments, the higher the substitution will score in the derived matrix.

Hopefully these scores are intupretable based on the biochemical properties of the amino acids. For example, R -> K scores 2, a  high score for a mismatch, but both amino acids are basic and so share a similar structure, presumably making a change more likely to occur and/or more likely to maintain the protein function and so not be harmful.

For your case, I don't have a good biological explanation for why E -> K scores positive, another author notes this is surprising, but ultimately the reason it is given a positive score is just because regardless of why, E -> K does apparently frequently occur, at least in the alignment data used to derive this scoring matrix. Incidentally, BLOSUM matrices are pretty ancient and based on few alignments by modern standards, but nonetheless appear to work quite well.


Need a alternative or more complex version of venn diagram in python for matching dna sequences

To get around limitations in using Venn diagrams for set overlaps, when there are more than three sets, back around 2013 or so I created something I called an Eulergrid plot (example at the bottom of the page), which an UpsetR plot appears to recapitulate, today.
The Eulergrid code I wrote was in a mix of Perl and R; the UpsetR plot code uses R. There appears to be a Python port of UpsetR that might work, as well.


See this post on the GATK forums.

It is normal to see values above 255.  The QUAL score is the posterior probability that all samples in your callset are homozygous reference. I quote Sheila's response from the thread posted above:

  The QUAL score is the Phred-scaled posterior of AC = 0. We use the AC priors and the PLs to get the likelihood of the data given each AC, then use those to get the posterior probability for each AC. From there, the calculation is 1 - Pr{AC > 0}.

Please see the GATK post in Methods and Algorithms for the nitty gritty details. Beware, Bayesian statistics abound.


Where can I find a GWAS data set about genotypes and phenotypes for Zebrafish and rice?

I think it will be more productive for you if you choose one of the public human datasets and practise your GWAS methods on that. A couple of examples examples:

1000 genomes -- use "population" as your trait of association. You'll find that you get association across the whole genome, so the analysis method would change from "what bits are associated" to "what are the bits that are most consistently associated"
NCI Genomic Data Commons -- open access (and controlled access) data associated with cancer


if I know the number of sequencing circles can I give this information to DESeq2?

You can include that info in the design, though you'd have to decide whether to include it as a number of just a factor.  I've never heard of anyone including PCR cycles.  I've heard of people including RIN scores, though I've never done it.
If the PCA looks alright, I wouldn't worry about.  I've got tons of examples of the lab people saying "I did a few more cycles on these few samples" and then you can't tell any difference in the PCA.  I think needing to do more cycles likely has more to do with variance in the prep than a real biological difference.


How does Cellminer's "Cross-correlations of transcripts, drugs, and microRNAs" work

In the article describing Cellminer says about the cross-correlations:

  ... a cross-correlation table of the resultant z-scores can be
  generated ...

On page 3503

And the z-scores are defined previously as:

  The tool output includes relative transcript intensity presented
  as z-scores...


Sort vcf by contig and position within contig

You could simply use bcftools sort for it:

$ bcftools sort input.vcf > output.vcf

If you really want to use bash only, you can do this:

$ grep "^#" input.vcf > output.vcf
$ grep -v "^#" input.vcf| sort -k1,1V -k2,2g >> output.vcf

The first command will write the header information to the new vcf file.
The second will sort by contig name and position and append the result to the new vcf file.

fin swimmer


How to upload and process a data set to Galaxy instance and retrieve results programmatically?

You best bet is bioblend, since it's a more convenient wrapper around the Galaxy API, but if you really want to you could use that directly. I should note that large files should probably still be uploaded via ftp or sftp rather than via the http(s) interface.


How to calculate RNA copies in qPCR using R?

Now that we've straightened out the data, it seems like the simplest solution is to just linearly interpolate the copy number from your standard curve.
When I read your standard curve (std) and your sample data (ct) into R, I see that the log10(Ct) fits the copy number fairly well, whereas the linear is a bit harder:

So we can fit a model to predict copy number based on Ct on the STD data:
> fit = lm(std$Log10.STD_copies. ~ std$STD_Ct)
> fit

Call:
lm(formula = std$Log10.STD_copies. ~ std$STD_Ct)

Coefficients:
(Intercept)   std$STD_Ct
    10.6847      -0.2357

# we plug those numbers into the linear equation,
# taking the mean of the two replicates for a best estimate.

> rowMeans(ct[,2:3]) * -.2357 + 10.6847
[1] 5.130429 5.210567 5.210567 4.595390 5.469837 6.582341
[7] 6.636553 6.391424

So those are log10(copy number) estimates, of a sort. You can get much more elaborate, of course.


Single-cell analysis to compare samples is a long a difficult process. There is very good documentation for 10x Genomics cellranger, the DropSeq Pipeline and the Seurat R package. These tools all have GitHub repositories and the authors are very responsive if you encounter issues.

Depending on the technology used to generate the data, you'll need to use either cellranger or DropSeq to process the FASTQ files. These are designed to account for the the different experimental design. Still, it is important to ensure that your samples have been processed in the same way to reduce batch effects.

This question is to open-ended to cover all of the details needed to perform the analysis but here is some of the main things to consider and the general process to perform such as analysis. You will need competence using the command-line (shell) and programming in R (or Python) to perform single-cell analysis. I recommend performing this analysis on a remote server as some steps are memory intensive.

The general overview to compare samples is:

Demultiplex the reads for different indicies. You need to do this even if you've only put one sample per lane (to remove the index from the reads). cellranger mkfastq or Illumina's bcl2fastq will do this.
Obtain a reference genome (FASTA) and gene annotation (GTF) for the species you are working with. You can prepare a reference transcriptome with cellranger mkgtf and cellranger mkref.
Run cellranger count or the DropSeq pipeline on each sample separately. These will both perform STAR (splice-aware) alignment of paired-end RNA-Seq reads and count UMIs for each cell-barcode and gene. Cells will be identified and filtered by an automatic UMI threshold (this can be changed by forcing with expected number of cells). This will return a gene-barcode matrix.
You need to aggregate results of cellranger runs for different samples with cellranger aggr. This will perform downsampling by default to normalise the number of reads between samples. It will also perform graph-based (Louvian) clustering on the combined samples and return HTML and Cloupe summary data to explore, in addition to a combined gene-barcode matrix for all samples.
Single-cell RNA-Seq experiments are subject to batch effects. The CCA method from Seurat and the MNN method from scran are both available in the respective R packages to account for this (with different approaches). Batch effect correction needs to be performed before downstream analysis to ensure comparisons are valid.

This methods require an overlap between the samples to estimate technical errors and account for them. You need to bear in mind that neither of these methods works unless:

some of the cells detected in each sample are expected to be the same cell type
biological replicates of each sample have been performed

Once you have performed this correction, then standard Seurat workflows such as those shown in the tutorials can be performed which will identify clusters. You can use these to identify clusters specific to each genotype. Bear in mind that some overlap should be expected and if there isn't one, batch effect correction will over-correct for this.


How are MACS2's narrow peak and broad peak algorithms different?

The key function is call_broadpeaks:

The description attached to the function says:

  This function try to find enriched regions within which, scores are
  continuously higher than a given cutoff for level 1, and link them
  using the gap above level 2 cutoff with a maximum length of
  lvl2_max_gap.

  scoring_function_s: symbols of functions to calculate
  score. 'p' for pscore, 'q' for qscore, 'f' for fold change, 's' for
  subtraction. for example: ['p', 'q']

  lvl1_cutoff_s:  list of cutoffs
  at highly enriched regions, corresponding to scoring functions.

  lvl2_cutoff_s:  list of cutoffs at less enriched regions,
  corresponding to scoring functions.

  min_length :minimum peak length,
  default 200.

  lvl1_max_gap   :  maximum gap to merge nearby enriched
  peaks, default 50.

  lvl2_max_gap   :  maximum length of linkage
  regions, default 400.

  Return both general PeakIO object for highly
  enriched regions and gapped broad regions in BroadPeakIO.

To give some basic explanation, the algorithm (briefly) appears to be as follows:

Two separate levels of peaks are called, level 1 (a higher pval, ie more significant) and level 2 (a lower pval). Level 1 is controlled by -p and level 2 is controlled by --broad-cutoff. When each peakset is called, they are immediately linked by the max gap parameter for each set.
Then, assuming that all level 1 peaks should be inside level 2 peaks (this is an explicit assumption by MACS2), the algorithm groups level 1 peaks inside level 2 peaks to output a broad peak.

...

This has a few implications:

The broad peak calls really come from the level 2 peaks alone (+ linking). The level 1 peak calls allow you to distinguish sub peaks (so that you can have gapped peaks).
Aside from the linking, the broad peak calls would be the same as narrow peak calls, if you called both with the same pval threshold (for example, if you set --broad-cutoff 0.1 in broad peak mode, and the -p 0.1 for narrow peak mode)


Installation of PRANK MSA in WLS Ubuntu 20.04 LTS

Their instructions are outdated and don't match what exists on their server. The problem you're running into is not at the cd step but even before that - at your first step even. Your wget fails because that file does not exist. Go to http://wasabiapp.org/download/ and see for yourself - there is no /prank/ directory in there.
I tried using curl, and it wrote the 404 page as a .tar.gz file, so if you did that, the tar step would error out.


Human Genome Alignment and Variant Calling Benchmarks

NA12878

My recommendation: download raw reads from Illumina BaseSpace and do alignment yourself. Aligning ~40X worth of human reads takes overnight with ~16 CPU cores.
You can acquire Platinum Genome data from ENA. There is a BAM available for download, but that was done by Illumina's own aligner against a wrong reference genome. You'd better do alignment by yourself. Another problem with this dataset is it is too good to be representative. Modern data are worse.
GIAB provides BAMs here:

ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/NIST_NA12878_HG001_HiSeq_300x/

It is not "typical", either.

SGDP

Instructions on downloading SGDP BAMs are available on this page, in particular this Word document. Warning: you need ~30TB disk space for these BAMs. I also heard it is tricky to download data.
There is apparently a mirror from cancer genomics cloud. I don't know how easy it is to access.


Oh you silly sausage!
This obviously means that you are running out of stack memory. Apparently...

freebayes is creating a single string containing the entire vcf header, which gets quite big when the reference has a lot of contigs with long names, causing stack space depletion.

Try increasing the maximum stack size (in this example to 128M) and running your command again.
ulimit -S -s 131072


I just received a reply from 1000Genomes regarding this. I'll post it in its entirety below:

  Looking at the example you mention, I find it difficult to come up with an
  interpretation of the information whereby the stated end seems to be correct,
  so believe that this may indeed be an error.

  Since the v4.0 was created, however, new versions of VCF have been introduced,
  improving and correcting the specification. The current version is v4.3
  (http://samtools.github.io/hts-specs/). I believe the first record shown on
  page 11 provides an accurate example of this type of deletion.

  I will update the web page to include this information.

So we can take this as official confirmation that we were all correct in suspecting the example was just wrong.


The following bit of python code should work:

#!/usr/bin/env python
import sys

lastTranscript = [None, None, None, []]  # ID, chrom, strand, [(start, end, score), ...]

def getID(s):
    """Parse out the ID attribute"""
    s = s.split(";")
    for k in s:
        if k.startswith("ID="):
            return k[3:]
    return None

def dumpLastTranscript():
    """Print the last transcript"""
    bounds = sorted(lastTranscript[3])
    print("{}\tgmapidx\tgene\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    print("{}\tgmapidx\ttranscript\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    for start, end, score in bounds:
        print("{}\tgmapidx\texon\t{}\t{}\t{}\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], start, end, score, lastTranscript[2], lastTranscript[0], lastTranscript[0]))

def handleLine(cols):
    """Handle a single line, appending the exon bounds to the previous if relevant"""
    ID = getID(cols[8])
    assert(ID is not None)
    if lastTranscript[0] is not None and lastTranscript[0] != ID:
        dumpLastTranscript()
        lastTranscript[3] = []
    lastTranscript[0] = ID
    lastTranscript[1] = cols[0]
    lastTranscript[2] = cols[6]
    lastTranscript[3].append((int(cols[3]), int(cols[4]), cols[5]))

f = open(sys.argv[1])
for line in f:
    if line.startswith("#"):
        continue
    cols = line.strip().split("\t")
    handleLine(cols)

dumpLastTranscript()
f.close()

If you saved that as gff2gtf.py and made it executable, the usage would be ./gff2gtf.py foo.gff3 > foo.gtf. With the example that you provided in your post, the result is:

2L  gmapidx gene    18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx transcript  18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx exon    18442664    18443024    79  -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
3R  gmapidx gene    15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx transcript  15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15853880    15855465    96  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15855529    15855742    97  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
X   gmapidx gene    14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx transcript  14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837470    14837753    92  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837810    14838142    93  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";

That at least looks like correct gtf 2.2 format to me.


How to perform bowtie2 analysis with slurm?

It seems you do not have bowtie installed or loaded. Use module load bowtie2 together when you load samtools.


dividing genome into non-overlapping windows using R

The question is a bit confusing to me, at it's core I understand you want non-overlapping windows across chromosomes.

One way to achieve is to use GenomicRanges::tileGenome function, which needs the chromosome lengths as input, e.g.

library(Biostrings)
mygenome <- readDNAStringSet(list.files(mypath,"mygenome.fa$",full=TRUE))
chrSizes <- width(mygenome)
names(chrSizes) <- names(mygenome)
print(chrSizes)
#     Chr1     Chr2     Chr3     Chr4     Chr5
# 18585056 19698289 23459830 26975502 30427671

This can than be fed into tileGenomes:

library(GenomicRanges)
bins   <- tileGenome(chrSizes, tilewidth=5e5, cut.last.tile.in.chrom=T)

This creates a GRange object which gives you a whole repertoire of post processing options, i.e. overlapping with SNPs or NGS reads (I use the latter for genome wide coverage plots):

print(bins)
# GRanges object with 240 ranges and 0 metadata columns:
#         seqnames               ranges strand
#            <Rle>            <IRanges>  <Rle>
#     [1]     Chr1   [      1,  500000]      *
#     [2]     Chr1   [ 500001, 1000000]      *
#     [3]     Chr1   [1000001, 1500000]      *
#     [4]     Chr1   [1500001, 2000000]      *
#     [5]     Chr1   [2000001, 2500000]      *
#     ...      ...                  ...    ...
#   [236]     Chr5 [28000001, 28500000]      *
#   [237]     Chr5 [28500001, 29000000]      *
#   [238]     Chr5 [29000001, 29500000]      *
#   [239]     Chr5 [29500001, 30000000]      *
#   [240]     Chr5 [30000001, 30427671]      *
#   -------
#   seqinfo: 5 sequences from an unspecified genome


Detecting whether genome is genetically engineered?

The definition of a GMO for example UK GMO regulatory bodies (HSE) is strictly the artificial introduction of foreign genetic material using targeted genetic engineering, such as splicing,  genome integration etc... Attenuation for example by serial passage the basis of the yellow fever vaccine, does not count as a GMO.

The simple answer is look at phylogenetic trees between the GMO strain against the total distribution of isolates. If large regions have been spliced the 'silent mutation rate' looks very different,.because the GMO parent is from a very different environment to the recipient.  The switch in replication rate affects the numbers of 'silent mutation', which is fast in RNA viruses.

If the GMO is a point mutation, ie limited amino acid change then I've actually done this with dengue virus and it is fairly easy because the novel mutation make the GMO 'stick out'.of the phylogeny against all natural mutation patterns. At this point the outlier is tracked back to Genbank where the relavant publication defines the engineered virus.

If you are inferring a conspiracy theory that is a different process and a different discussion, particularly if it is Ebolavirus which is a Category 4 pathogen. The investigation is again within the phylogeny, but in this case the answer is always the isolate is a natural zoonose from a sylvatic (wild) transmission cycle which has spilled over into humans.


Improve scRNA-seq dataset for further analysis

I don't think you can conclude that the dataset is terrible based on that PCA. Depending on the specific protocol, each scRNA-seq dataset is going to be very different. Unlike bulk RNA-seq where all the samples are going to be of very similar quality, individual cells will be highly variable. For example, the most basic QC metric is the number of reads. If some of your bulk RNA-seq samples have 10X more reads than other samples, that's a cause for concern. In scRNA-seq experiments, that's not at all surprising. The PCA exposes the major sources of variation. In scRNA-seq experiment, those are probably going to be technical rather than biological, so the output may not make much sense. Furthermore, even if all the cells are of perfect quality, PCA is not really optimal for sparse data with many dropouts like you see with scRNA-seq.

Regardless of PCA, those library sizes seem very small. Based on your plot, there are barely any cells with over 10K reads. I understand this is C elegans, but that seems somewhat low. That may be why you are having trouble with normalization. It would help with understanding the data if you changed the X-axis on those plots so that you can actually see where the peak is. Right now, it looks like 0, but that can't be the case.


How can I extract information from .sam files?

Use samtools flagstat with option -O tsv:

Using -O tsv selects a tab-separated values format that can easily be imported into spreadsheet software. In this format the first column contains the values for QC-passed reads, the second column has the values for QC-failed reads and the third contains the category names.

SEE ALSO:
What Does Samtools Flagstat Results Mean?
INSTALLATION:
Some old versions of samtools do not have option -O tsv. To install a recent versions of samtools, use either:
conda (macOS, Linux) (but see also: Tutorial: installing samtools with conda - overcoming two common errors):
conda create --name samtools samtools

or brew (macOS):
brew install samtools


Build a Pubmed query given long gene list


How to identify latent variables in single-cell RNA-Seq data

I think to identify latent variables, PCA is probably not going to work. NMF might be worth trying. You might want to check out a method called consensus NMF (cNMF) (https://elifesciences.org/articles/43803) published in eLife. The author described some latent variables corresponds to cell types, while others corresponds to cell-cycle or metabolic states of cells.


Counting repeated kmers sequences that match at least x % of reads sequence

The following javascript does what you want. You need node.js to run it. It should be easy to translate the code to Python. I have not carefully tested it. Use with caution.

EDIT (response to new comments): the program has been changed to compute the sum of lengths. Note that only a tandem repeat of length k*2 or longer is counted. For example in sequence ATTGGGATTGGGATTcGGGATTGGG, the script returns 15, because the second is not long enough. Counting length smaller than k*2 will be slower and more complex, but is often not the right thing to do. For example, in ATTGGGATTGGGcccccccgaaatcgatagcatcgaGGGATTGcgatc, we would not count the second GGGATTG as a tandem repeat.

On performance, for an input string of length $l$, the time complexity of this script is $O(l)$ associated with a small constant. It will be efficient and adequate to practical uses.

If you use TRF, you should run it on your sequences and then test if any tandem repeats it finds match GGGAAT – the procedure will be similar to my script, except that TRF will also find impure repeats.

// given the repeat unit length _k_, find all tandem repeats of 2*k or longer
function trf_k(k, str)
{
    var streak = 0, a = [];
    for (var i = k; i <= str.length; ++i) {
        if (i < str.length && str[i] == str[i - k]) {
            ++streak;
        } else {
            if (streak >= k)
                a.push([i - streak - k, i]);
            streak = 0;
        }
    }
    return a;
}
// find the sum of tandem repeat length with _kmer_ being the repeat unit
function trf_kmer(kmer, str)
{
    var a = trf_k(kmer.length, str), sum = 0;
    for (var i = 0; i < a.length; ++i)
        if (str.substring(a[i][0], a[i][1]).indexOf(kmer) >= 0)
            sum += a[i][1] - a[i][0];
    return sum;
}

var str = 'xATTGGGATTGGGATTGxGATTGGGATTGGGATTGGGATTGxx';
console.log(trf_kmer('GGGATT', str));


How to QC overamplified shotgun library?

FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) is a good start: it will show you observed/expected ratios for the most overrepresented sequences, and running it on a few datasets you know to be high quality should give you intuition for what looks good and what doesn't.

My labmate also produced something called preseq (http://smithlabresearch.org/software/preseq/) that uses the library complexity of an initial sequencing experiment to predict how many more unique reads you would get with additional sequencing of the same library. It could give you some intuition for whether or not it is worth sequencing your library deeper to get the information you need.


Software recommendation: find DNA sequence distribution over entire transcript

If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio:
library(EnsDb.Mmusculus.v79)
library(ggbio)
library(biomaRt)
library(stringr)
library(dplyr)

# plot gene model with ggbio
gene_model <- autoplot(EnsDb.Mmusculus.v79, ~ symbol == "Grin1")

ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

# You need the full length gene sequence to match the gene model
seq <- biomaRt::getSequence(id = "Grin1",
                         type = "mgi_symbol",
                         seqType = "gene_exon_intron",
                         mart = ensembl)

# retrieve position of gene in genome to match up with gene model
gene_start=biomaRt::getBM(attributes=c("start_position"), filters = 'external_gene_name', values = 'Grin1', mart=ensembl)

# use stringR to find positions of exact matches
motif <- "ggcc"
motif_positions <- as.data.frame(str_locate_all(seq$gene_exon_intron, toupper(motif))) %>%
  mutate(start=start+gene_start[[1]]-1)

# make plot similar to Llopis's answer
motif_plot <- ggplot(motif_positions, aes(x=start)) +
    geom_point(y=1) +
    scale_y_continuous(limits=c(0,1))

# combine plots
tracks(motif= motifs, gene_model=gene_model)

This gives all motif matches for the gene sequence, not just the transcripts. It would be quite straightforward to filter out positions that overlap introns


Is there a JSON-based genomic feature format?

JVCF is a JSON schema for describing genetic variants. BioJSON is a JSON schema for multiple sequence alignments. The output of mygene.info's gene annotation web service is JSON-formatted.

You might build one of your own schemas from any or all of these, as inspiration for structure and content.

Consider that most genomics formats are historically tabular, because they are not hierarchical, because they originated prior to the advent of JavaScript as a common programming language (and subsequent use of JSON to serialize data), and because a lot of stream-oriented UNIX tools read in data and process it in ways that are natural to tabular structure.

Some tabular results can be hierarchical, but use inefficiencies in field values to apply hierarchy. JSON is inefficient for storing tabular data, but it is absolutely wonderful for hierarchical data. This difference is perhaps worth considering.

One could think of storing a gene as a master object containing various haplotypes, exons, introns, features, associated variants, associated TF binding sites, gene ontologies, etc. Ordering or prioritization of such items would be possible with the use of sorted lists within a gene object. The format is as extensible as it is a matter of adding a property, list, or object to the schema.

To gain adoption or "buy-in", build tools which publish it to users, and tools that consume what users get in order to do things, along with conversion and validation utilities.


Inflate operation failed: progress temporarily not possible, or in() / out() returned an error

My best guess would be a corruption during download as:
samtools quickcheck -qvvvv ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
verbosity set to 4
checking ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
opened ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam is sequence data
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam has 196 targets in header.
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam has good EOF block.

Plus I can successfully run idxstats from remote and fetch individual chromosomes without error via samtools view. Try redownloading before digging deeper.


Based on the way you've called the command (using -xl-sn), I assume you are using GATK 4 and not GATK 3.

If this is the case, the documentation for the -xl-sn flag is available here: https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_variantutils_SelectVariants.php#--exclude-sample-name.

By default, the -xl-sn SAMPLE means "exclude the sample called SAMPLE from my VCF". In your case, you're providing an entire list of samples, in a file, which is not the default behaviour. In order to tell GATK that you're providing a file and not a sample name, you need to ensure the file you provide has the extension .args.

From the documentation:

  File names must use the extension ".args"


Demultiplexing and preprocessing for custom single-nucleus Drop-seq data

I didn't notice the sentence:

  Paired-end sequencing reads were processed largely as described (http://mccarrolllab.com/wp-content/uploads/2016/03/Drop-seqAlignmentCookbookv1.2Jan2016.pdf ), with additional correction steps.

This should provide instructions on how to lead the preprocessing and mapping of data.


How to resolve in snakemake error : "Target rules may not contain wildcards."

The problem is that you need a master rule that requires all of your desired outputs as inputs, in your case it would be :

rule all:
    input:
        expand("data/{sp}/genome.fa.gz", sp=species.split(' '))

You'll also need separate download link inputs for each species. You could make a separate download_table.tsv for each species, but it would probably be easier to make a config file with this information, and add a params keyword to your rule. Something like:

rule download:
   params:
       url=config['locations']['sp']
   output :
       "data/{sp}/genome.fa.gz"
   shell :
       "scripts/download.sh {params.url}"


How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq?

Actually this is one of the main problems you have when analyzing scRNA-seq data, and there is no established method for dealing with this. Different (dedicated) algorithms deal with it in different ways, but mostly you rely on how good the error modelling of your software is (a great read is the review by Wagner, Regev & Yosef, esp. the section on "False negatives and overamplification"). There are a couple of options:

You can impute values, i.e. fill in the gaps on technical zeros. CIDR and scImpute do it directly. MAGIC and ZIFA project cells into a lower-dimensional space and use their similarity there to decide how to fill in the blanks.
Some people straight up exclude genes that are expressed in very low numbers. I can't give you citations off the top of my head, but many trajectory inference algorithms like monocle2 and SLICER have heuristics to choose informative genes for their analysis.
If the method you use for analysis doesn't model gene expression explicitly but uses some other distance method to quantify similarity between cells (like cosine distance, euclidean distance, correlation), then the noise introduced by dropout can be covered by the signal of genes that are highly expressed. Note that this is dangerous, as genes that are highly expressed are not necessarily informative.
ERCC spike ins can help you reduce technical noise, but I am not familiar with the Chromium protocol so maybe it doesn't apply there (?)

since we are speaking about noise, you might consider using a protocol with unique molecular identifiers. They remove the amplification errors almost completely, at least for the transcripts that you capture...

EDIT: Also, I would highly recommend using something more advanced than PCA to do the analysis. Software like the above-mentioned Monocle or destiny is easy to operate and increases the power of your analysis considerably.


Humans are diploid, so you can expect to see up to 2*N (2*68=136) alleles, so N_DATA is the number of observances of that allele.


RNA_Seq Analysis in R, propmapped function issue

There must be a problem in the BAM files you are using. Here are steps I would do to figure where exactly is the problem.

Check the data you have downloaded. The tutorial is using 12 fastq.gz files (SRR1552444 - SRR1552455). You can check on sra how big they are supposed to be (For example SRR1552444 is supposed to be an ~1.7Gb file) and if your downloaded files are about the same. Check if the files are fine (In the terminal you can type zcat data/SRR1552444.fastq.gz | head to check the first few reads). If the files are 0 sized or somehow corrupted, you need to redownload them.

Check if your mapping reference looks alright (size wise, look at the first few lines) and if the reference is indexed (there is a .fai in the directory).

Check the size of the BAM files, they also need to be in the range of a few hundred Mb to a few Gb. I expect these will be corrupted for sure, otherwise, you would get your numbers. You might also want to check if they are indexed (although this should not be a problem, I suspect they will get indexed automatically). If there is nothing wrong the bamfiles, you might want to check another tool to summarise a bamfile (like samtools flagstat) to check if at the end the problem really is in propmapped. If that's the case.

This should at least help you to zoom to the core of your problem if not resolve it.


how do I convert a vcf file for a GWAS study (using R package vcfR)

If you want to do GWAS or mendelian randomisation, you can do it with Plink (v2) which should be faster than R.


Y Chromosome Aligned Reads in scATAC-seq data from a female-derived cell line?

There are homologous regions between X an Y chromosomes: https://en.wikipedia.org/wiki/Pseudoautosomal_region

It is therefore normal to have some female-derived reads mapping in Y chromosome.

You should probably check what proportion of such reads fall in other parts of the Y chromosome than pseudoautosomal regions.


How to trim adapter sequences from GSE65360 in order to map the reads?

First the Nextera adapters and the custom barcode adapters overlap each other

Nextera 1                                     TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
Barcode 1   GATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT

Nextera 2                                    GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
Barcode 2    CAAGCAGAAGACGGCATACGAGATACCCAGCAGTCTCGTGGGCTCGGAGATGTG

Looking at some read 1s and read 2s you were correct that the adapters appear at the 3' end of the read reverse complimented (The |> points to the adapter in the sequence).  You can see that CTGTCTCTTATACACA would be the ideal sequence to use for the 3' end of the reads.  Since the Nextera adapter and the barcode overlap you only need to use the one that is adjacent to the insert.

SRR1780164_2.fastq: GCCCGAGGGGACGCTGCCCCACGTGCTCCAGG |> CTGTCTCTTATACACATCTGACGCTGCCGACGATGCACGAAGTG
SRR1780164_1.fastq: GGGCTGGACAGCGGCGGGGACGCGGACCTGCAG |> CTGTCTCTTATACACATCTCCGAGCCCACGAGACACAAAGTGA

Of the 64k reads there are ~11k with the Nextera adapter at the 3' end based on a rough sequence grep

On the 5' end of the read there is some Nextera adapter present but it is not reverse complimented and there are only a handful of reads with the adapter present ~250 total between read 1 and read 2.

SRR1780164_2.fastq: TGGGCTCGGAGATGTGTATAAGAGACAG <| CTCCCCCGCCTCCAGCATCCGGGCGAGGTAGTGCATCGACGCGTCCAC

I used this command:

cutadapt --minimum-length=20 -A CTGTCTCTTATACACA -a CTGTCTCTTATACACA -G GATGTGTATAAGAGACAG
-g GATGTGTATAAGAGACAG -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz
SRR1780164_1.fastq SRR1780164_2.fastq

Which gives a trimming result of

   === Summary ===

Total read pairs processed:             64,832
  Read 1 with adapter:                  16,721 (25.8%)
  Read 2 with adapter:                  16,996 (26.2%)
Pairs that were too short:                  46 (0.1%)
Pairs written (passing filters):        64,786 (99.9%)

I think this best represents the library construct and the trimming of the adapter sequences


Are duplicate variants against the VCF standard?

No. The only point the specification makes is that:

  It is permitted to have multiple records with the same POS

I suspect this is intended to let the file list different variants occurring at the same position, but this incidentally allows duplicate variants.

This said, you probably don't want duplicate variants, because it will mess with any statistics you apply over the variants, and will unnecessarily clutter the file and slow down processing, if this happens a lot. You might also get variant annotations on one variant but not another, which can become a problem depending on how you query your variants.


From your example we are comparing samples not genes, I think we need to run the test row wise, something like this:

pvalueGenes <- sapply(colnames(biomarker),
                      function(i) t.test(biomarker[i, ],
                                         immune[i, ], paired = TRUE)$p.value)

names(pvalueGenes)[ pvalueGenes < 0.05 ]


Is there a more elegant solution to the bwa-mem: paired reads have different names error?

Assuming the read name lines only contain things like @JQ714243-920-1 and nothing else, then sed -i '1~4 s/-[12]$//' file.fastq (N.B., example updated a bit for safety, thanks to @terdon) will remove the -1 and -2 suffixes.

However, you should try to figure out why the read names were munged into that form to begin with (at least Illumina machines won't produce that, I don't know about others).


use same output in two processes in nextflow dsl2

I have solved this problem in the past by duplicating the output channel of the desired process in order to feed the same output into multiple downstream processes.


To set the ball rolling, I tried to implement it manually with your files. Overall, there are three steps:

Locate NITNLC (or HAPATV) in protein P0DTC2.
Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1.
Print the amino acid and the DNA codon of both proteins from step 2.

It works, but only for the first 60 amino acids. I wonder if it's due to the amino sequences in AMINOOO_seq_removed.fasta repeat itself every 60 acids. But why?
#For example, the first three lines of protein P0DTC2
>P0DTC2
MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS
NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV
NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE
...

Step 0: Read the files.
library('seqinr')

#align file containing protein sequences
count_added <- read.alignment('count_added_.clustal_num', format='clustal')
names(count_added$seq) <- count_added$nam

#DNA sequences
rev3.aln <- read.alignment('rev3.aln', format='fasta')
names(rev3.aln$seq) <- rev3.aln$nam

Step 1: Locate NITNLC (or HAPATV) in protein P0DTC2. Again, there are two repeats of NITNLC 60 acids apart (at 807 and 870).
locate <- function(seq, find)
  {address <- gregexpr(paste(strsplit(find, '')[[1]], collapse='[^a-z]*'), seq)
  #substr(seq, address[[1]][1], address[[1]][1]+attr(address[[1]], 'match.length')[1]-1)

  return(list(start=as.numeric(address[[1]]),
             end=as.numeric(address[[1]] + attr(address[[1]], 'match.length') - 1)))
  }

locate(seq=count_added$seq[['P0DTC2']], find='nitnlc')
#start
#807 870
#end
#812 875
locate(seq=count_added$seq[['P0DTC2']], find='hapatv')
#start
#1212 1279
#end
#1217 1285

Step 2: Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1. Instead of using the range 812 - 1279, I chose range 1 - 20 for demonstration purposes.
compare <- function(seq1, seq2, after=0, before=100000)
  {seq1_ = strsplit(seq1, '')[[1]]
  seq2_ = strsplit(seq2, '')[[1]]
  ind = which(seq1_ != seq2_ & grepl('[a-z]',seq1_) & grepl('[a-z]',seq2_))
  ind = ind[ind>after & ind<before]
  #seq1_[ind[1]]
  #seq2_[ind[1]]

  return(ind)
  }
compare(seq1=count_added$seq[['P0DTC2']], seq2=count_added$seq[['K9N5Q8']], after=1, before=20)
# [1] 5  7  8  9 10 13 14

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34

Step 3: Print the amino acid and the DNA codon of both proteins from step 2. Note that ind is based on the align file.
print_amino_codon <- function(ind, seq, seq_gene)
  {locate_amino <- gregexpr('[a-z]', seq)[[1]]
  if (!ind %in% locate_amino) return(NA)
  ind2 = match(ind, locate_amino)

  return(c(amino=substr(seq, ind, ind), codon=substr(seq_gene, ind2*3-2, ind2*3)))
  }

codon(ind=5, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "v" "atg"
codon(ind=6, vseq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "f" "ttc"
codon(ind=7, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "l" "ttg"

codon(ind=5, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "m" "atg"
codon(ind=6, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "f" "ttt"
codon(ind=7, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "v" "gtt"

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34
#<print>         ^^^

#K9N5Q8 gene codon
#gtg ttt cta ctg atg ttc ttg tta aca
#                ^^^ ^^^ ^^^
#P0DTC2 gene codon
#atg ttt gtt ttt ctt
#^^^ ^^^ ^^^


Check your webserver functions from the top down.  You can do this from any WEB browser, even one that is on the webserver system itself.  For example, does this link function:
  http://webserver-id/
and provide you with a home page for that WEB service ?  If not, then there
is not apache WEB server on that system, or is mis-configured.

Next, depending upon how that Apache is configured, try this URL
  http://webserver-id/~username/
which would be a normal way to get to user directories via Apache.
This would actually be looking for a special directory in your HOME
directory, perhaps something like $HOME/public_html/
You would need to place something in public_html/ to access, for example
a tiny WEB page:

<title>my home directory</title>

hello world

If that works, then you need to place your hub files in this
public_html/ directory.


Any way to quantify the variation of genes that expressed in Affymetrix expression data?

You can use the following code to calculate the coefficient of variation:

# expr is your expression matrix.
SD <- apply(expr, 1, sd)
CV <- sqrt(exp(SD^2) - 1)

It might be implemented in some package but it is so brief that you can write it again yourself. Then you can filter out those that are below certain percentage of the distribution of CV (like 0.1)

expr[CV > quantile(CV, probs = 0.1), ]


How to loop multiple function in shell script?

You don't need a loop. You can do the whole thing with a simple awk one-liner:

awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta

I ran this on your files and got:

$ for file in gene*; do echo "=== File: $file ==="; cat $file;  done
=== File: gene1.fasta ===
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT
=== File: gene2.fasta ===
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA
=== File: gene3.fasta ===
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Explanation

-vRS='>' : set the record separator to >. This will make awk consider every fasta sequence as a single "line".
FNR>1{foo} : FNR is the line number of the current file. Because > is the first character of the entry, the first "line" will always be empty, so it needs t obe skipped.
printf ">%s",$0 > "gene"FNR-1".fasta" : this does the work. It prints a > followed by the current "line" ($0, this is the entire sequence + the header) and redirects it into a file called gene, then the current line number -1 (-1 is needed because of the issue mentioned above) and a .fasta extension.

=====

Note that org*fasta will be expanded with the default collation order. This means that org10.fasta will be processed before org1.fasta:

$ printf "org1.fasta\norg10.fasta\n" | sort
org10.fasta
org1.fasta

To avoid this and have the files in numerical order, you can use brace expansion instead of the *. So, if you know your files are all orgN.fasta where N can take values from 1 to 20, do this:

awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org{1..20}.fasta


how to write the script to know each sample's cluster number?

Assuming you are following the basic seurat guide:

pbmc@ident

You can access the 'ident' slot of your Seurat object and you will get a vector of the cluster IDs and the names will be the respective cells.

head(seuratobj@ident)
Cell1A7 Cell1C4 Cell1D1 Cell1D2 Cell1D3 Cell1D5
0 3 3 1 3 1
Levels: 0 1 2 3 4


Is there a standard definition for "assembly polishing"?

In long-read assembly, "polish" refers to the step to improve the base accuracy of contig sequences. I believe the terminology was originated from the HGAP paper:

  The final consensus–calling algorithm Quiver, which takes into account all of the underlying data and the raw quality values inherent to SMRT sequencing, then polishes the assembly for final consensus accuracies in excess of 99.999% (QV of >50)...

Pilon was published more than a year later and in fact, neither its paper nor README mentions word "polish".

By the way, "error correction" refers to the step to fix sequencing errors in raw reads, not in contigs. In Canu/HGAP/Falcon, that step is done by pbdagcon or falcon_sense.


Techniques for analyzing and quantifying sample bleed through in genotyping with Illumina

Are you using dual-index barcodes with different barcodes at each end? There is a known phenomena of "index switching" that occurs in Illumina reads. One way to control for this is to add in unamplified water samples that have unused barcode combinations where each of the indexes is shared with another sample.

More information can be found on this SeqAnswers thread. Here are some useful associated links:

Signal-spreading preprint -- detailed investigation of the phenomenom, including methods, graphs, and associated data
Crossblock -- tool by Brian Bushnell that attempts to remove cross contamination
Summary by James Hadfield -- mentions a few fixes suggested by Illumina, many centred around not using barcode combinations that could conflict
Illumina's quick-fire response -- they've known about it for a long time, and point out that it shouldn't be a problem for most purposes

I have suspicions that at least some of the index switching is being caused by ligation of two separate barcoded fragments during the sample prep (which definitely happens during sample prep for nanopore), but I haven't done any experimentation to work out if this is also the case for Illumina (due to money/time/goal constraints). The experimental design would be to take a prepared Illumina library just prior to sequencing, ligate on ONT adapters, then sequence on both the ONT and Illumina machines to see if chimeric reads discovered via nanopore sufficiently matches the index switching observed via short-read sequencing.


sum of products of mapq and mapped bases for each read in a from a BAM file

BAM files are stored line by line, so your only option is to process it that way. In python with pysam:

import pysam

bam = pysam.AlignmentFile("something.bam")
tot = 0
for read in bam.fetch():
    if read.is_unmapped:
        continue
    tl = 0
    for op, l in read.cigartuples:
        if op in [0, 7, 8]:
            tl += l
    tot += tl * read.mapping_quality
bam.close()
print("the total is {}".format(tot))


DESeq2 multiple treatments, multiple time points, multiple cell lines

To compare a subset of samples to another subset of samples, make one column of column data which will distinguish them all the way you want (so probably join cell line, treatment, timepoint), and make that column your design, and then when you call results, you specify the pair you want to contrast, one pair at a time.


match two files and check if position of chrom is within exon region

Despite an awk solution is favored by the OP, here is an R solution. (The task, especially dealing with exon coordinates in file2, is rather complex for a straightforward awk solution, at least for me).

res <- data.frame() # empty df for storing results

for (i in 1:nrow(file1)) {

  chr <- file1[i,1] # get chr
  pos <- file1[i,2] # get position of interest
  print("pos")
  print(pos)

  if(chr %in% file2$chrom) {

    file2_subset <- file2[file2$chrom == chr,] # subset file2

    for(j in 1:nrow(file2_subset)) {

      # create exon start - exon end pairs

      exon_start <- file2_subset[j,2]
      exon_start <- as.numeric(unlist(strsplit(exon_start,",")))
      print(exon_start)

      exon_end <- file2_subset[j,3]
      exon_end <- as.numeric(unlist(strsplit(exon_end,",")))
      print(exon_end)

      # if position of interest lie within exons
      # need to be careful if start/end position are counted as exon

      if(any(pos >= exon_start & pos <= exon_end)) {
        print(file1[i,])
        res <- rbind(res, file1[i,]) # write line i from file1 if satisfies above conds
      }
    }
  }
}

> res
  CHROM    POS REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2          INFO
6  chr1 790590   C   T    0/1    0/1       0|0       0/1       0/1  0|0  1|1 AN=2184;AC=84


The counts are "reads" for single-end datasets and "fragments" for paired-end datasets. In other words, they're what featureCounts and htseq-count produce. The "normalized counts" that you'll be able to access in DESeq2 are, as aptly named, simply normalized versions such that values are comparable across samples.


State of the art in predicting Translation Initiation Sites

Ribosome Profiling experiments profile the positions of ribosome-protected mRNA. a cycloheximide treatment fixes the translating ribosomes. A naive approach to define Translation Initiation Sites (TIS) would be to look at positions where you see distinct 'peaks' in the distribution of these fragments since the ribosome density will be higher at the start of the CDS (at the TIS).

TISDB is a database comprised of GTI-Seq global mapping of TIS codons at nearly single-nucleotide resolution. GTI-Seq uses two 'fixers' that fix ribosomes where one of them specifically fixes the initiating ribosomes.


Concatenate multiple fasta files into one file per unique entry

You can do this simply enough in awk:

$ awk -F'_' '{
               split(FILENAME,n,/\./);
               if(/^>/){
                  sample=$1;
                  sub(/>/,"",sample);
                  $0=">"n[1]
               }
               print > sample".fasta"
             }' gene1.fasta gene2.fasta ... geneN.fasta

If I run this on the two files you gave, I get:

$ cat sample1.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

$ cat sample2.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

Explanation

-F'_' : set the field separator to _ as an easy hack that lets us use $1 for the sample name (sample_NNN, so sample will be $1 and NNN will be $2).
split(FILENAME,n,/\./); : FILENAME always holds the name of the file currently being processed. Here, we are splitting it on . to get the gene's name. This assumes all of your files are named $geneName.fasta. The split creates the array n whose 1st element, n[1] is the gene name.
if(/^>/){ : if this is a fasta header line.
sample=$1 : set the variable sample to the sample's name. See the first note about -F'_'.
sub(/>/,"",sample); : remove the > from the sample variable. We don't want that in the file name.
$0=">"n[1] : set the current line to be > and the gene name.
print > sample".fasta" : print each (modified) line into the appropriate file name.


Clustering of gene co-expression network by igraph R package

For the first part, do you mean that the file is too large to be run on your computer?

For the second, if I understood correctly, you can use igraph or ggnet2 to color code the nodes based on the conditions you have and see if communities are formed based on condition. This would be a simple way to do this kind of thing.

I personally prefer using ggnet2 as I found it more user friendly. What you can do if you have an adjacency matrix, say adj_mat is to create a color vector based on the conditions for each node in the network, where length(color_vec) = ncol(adj_mat) = nrow(adj_mat) - since it is a square matrix. Following that, you could do something like this:

library(ggplot2)
library(network)

neto <- network(adj_mat, directed = FALSE, loops = FALSE, hyper = FALSE)
#assuming that it is undirected here and a simple network diagram.
#run ?network for more details in case your analysis is to be more complex

color_vec <- as.factor(color_vec)

col_pal <- c("blue", "red", "green") #assuming you only have 3 colors..

names(col_pal) <- levels(color_vec)
neto %v% "color" <- as.character(col_vec)

ggnet2(neto, node.size = .9, color = "color", palette = col_pal,
       edge.size = .05, edge.color = "white", mode = "adj")


Can anyone point me towards a HiC Dataset with domains of high confidence?

I guess your best chance is one of the in-situ HiC papers from the Lieberman-Aiden lab. There are not too many HiC datasets with 10kb-or-better resolution out there for human. I would check this one:
A 3D Map of the Human Genome at Kilobase Resolution Reveals Principles of Chromatin Looping


There's a faster way to download 1000 Genomes phase 3?

You are using the command for older versions. In the page you link to, they give two commands, one for the newer and one for the older versions of aspera. You are using the one for the older. The only difference is the keyfile, but what you want is:

ascp -i /home/ibmuser/.aspera/connect/etc/asperaweb_id_dsa.openssh -Tr -Q -l 100M \
 -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi ./

At least, I just tried this a second ago with the aspera bundle I downloaded today and that's the command that worked: asperaweb_id_dsa.openssh and not asperaweb_id_dsa.putty.


DoHeatmap Hierarchical Clustering Seurat

You need to order the marker matrix (e.g. by avg_logFC) before calling DoHeatMap.

library(dplyr)
all.markers <- FindAllMarkers(object = obj)
top20 <- all.markers %>% group_by(cluster) %>% top_n(20, avg_logFC)
DoHeatmap(object = obj, genes.use = top20$gene, slim.col.label = TRUE, remove.key = TRUE)


How can I extract only insertions from a VCF file?

A one-liner:

zcat my.vcf.gz |
  perl -ane '$x=0;for $y (split(",",$F[4])){$x=1 if length($y)>length($F[3])}print if /^#/||$x'

or equivalently

zcat my.vcf.gz |
  perl -ane '$x=0;map{$x=1 if length>length($F[3])}split(",",$F[4]);print if /^#/||$x'

For simple VCF operations, I generally recommend to write a script. This may be faster than those using heavy libraries. With a script, you only parse fields you care about; most libraries unnecessarily parse every field.

On a related note, I recommend not to decompose multi-allelic sites unless necessary. Decomposing is tricky, makes VCF harder to parse and to understand and may be a potential source of errors. Here is an example:

#CHROM POS ID REF  ALT               QUAL FILTER  INFO  FORMAT  S1  S2  S3  S4
11     101 .  GCGT G,GCGA,GTGA,CCGT  199  PASS    .     GT      0/1 1/2 2/3 2/4

vt decompose+normalize produces the following VCF:

#CHROM POS ID REF  ALT QUAL FILTER INFO FORMAT  S1   S2   S3   S4
11     101 .  GCGT G   199  PASS   .    GT      0/1  1/.  ./.  ./.
11     101 .  G    C   199  PASS   .    GT      0/.  ./.  ./.  ./1
11     102 .  CGT  TGA 199  PASS   .    GT      0/.  ./.  ./1  ./.
11     104 .  T    A   199  PASS   .    GT      0/.  ./1  1/.  1/.

In theory, you can reconstruct the original VCF from this output. However, it is very challenging for a program to do that. When you compute allele frequency line-by-line, this VCF will give you wrong results. bcftools norm -m- replaces "." with "0". You can get a correct ALT allele frequency from the bcftools output, but a wrong REF allele frequency. Furthermore, vt is also imperfect in that "CGT=>TGA" is not decomposed.

My preferred output is:

#CHROM POS    ID     REF    ALT    QUAL   FILTER INFO   FORMAT S1     S2     S3     S4
11     101    .      GCGT   G,<M>  0      .      .      GT     0/1    1/2    2/2    2/2
11     101    .      G      C,<M>  0      .      .      GT     0/2    2/0    0/0    0/1
11     102    .      C      T,<M>  0      .      .      GT     0/2    2/0    0/1    0/0
11     104    .      T      A,<M>  0      .      .      GT     0/2    2/1    1/1    1/0

Here we use a symbolic allele <M> to represent "another ALT allele". You can calculate the allele frequency by looking at one line, and won't confuse other ALT alleles with REF. bgt can produce such a VCF indirectly. However, it discards all INFO, so is not a practical solution, either.

In summary, it is very difficult to decompose multi-allelic sites. When you get decomposing wrong, your downstream analyses may be inaccurate. Decomposition should be used with caution.


Get regions' information from DNA sequence data (bsgenome.hsapiens.ucsc.hg19)

You can get specific sequences efficiently using the following packages:

library(GenomicFeatures)
library(Biostrings)
library(tidyverse)
library(data.table)
library(dtplyr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

Methylation_data <- tribble(
  ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
  "chr1",15864,15866,FALSE,894,"+",
  "chr1",534241,534243,FALSE,921,"-",
  "chr1",710096,710098,FALSE,729,"+") # Arrange your dataset

Methylation_GRanges <- GRanges(  # Construct a GRanges objects
  seqnames = Methylation_data %>% pull(seqnames), # First pull in the chromosome locations
  IRanges(start = (Methylation_data %>% pull(start))-50, # then the start position minus 50 nt
          end = (Methylation_data %>% pull(end))+50), # the end position plus 50 nt
  strand = Methylation_data %>% pull(strand) # the strand your sequence is located on
)
> getSeq(Hsapiens,Methylation_GRanges)
A DNAStringSet instance of length 3
width seq
[1]   103 GCTGCTGCTTCTCCAGCTTTCGCTCCTTCATGCTGCGCAGC...CTTGGCGGATGGACTCTAGCAGAGTGGCCAGCCACCGGAG
[2]   103 ACAGTTCCAATGTAATCAGAGAGAACATCACACACACACCA...ACCTGAGCAGCACTCTGCAAAGCTGTCAAGGCGGTGAAGC
[3]   103 GATAAAATAAAGCTTAGATTGGAAAAAATATTTAAGATTCT...GGAAGCTGAGTAATTGTATGTTCAAATACTTGCAAAACAT

To write the sequences to a file, use writeXStringSet()

For the second part of what you need, GenomicFeatures has a number of functions that obtain the annotated sequences of the whole genome by its type:

transcriptsBy(x, by=c("gene", "exon", "cds"), ...)
exonsBy(x, by=c("tx", "gene"), ...)
cdsBy(x, by=c("tx", "gene"), ...)
intronsByTranscript(x, ...)
fiveUTRsByTranscript(x, ...)
threeUTRsByTranscript(x, ...)

And this is how you use one of them to find whether your sequences overlap a particular type of DNA:

Find_overlap <- function(Ref_Data,
                         Overlap_Message){

  Ref_Data <- as.data.table(Ref_Data) # Convert to data.table for efficiency

  Methylation_data <- tribble(
    ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
    "chr1",15864,15866,FALSE,894,"+",
    "chr1",534241,534243,FALSE,921,"-",
    "chr1",710096,710098,FALSE,729,"+") %>%
    mutate(Overlap = "-",
           start = start-50,
           end = end+50) %>%
    as.data.table()

  Methylation_data[Ref_Data,
                   Overlap := Overlap_Message,
                   on = .(start >= start,
                          end <=  end)] # wherever the range overlaps, assign an identifier to it
}

Methylated_in_5UTR <- fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene) %>%
  Find_overlap(., "5' UTR") %>%
  print

seqnames  start    end Logical Methylation Level strand Overlap
1:     chr1  15814  15916   FALSE               894      +       -
2:     chr1 534191 534293   FALSE               921      -       -
3:     chr1 710046 710148   FALSE               729      +  5' UTR
```


How can I speed up INDEL calling/correction on BAM files?

I wasn't aware of the samtools subsampling when I had this problem a couple of years ago, so ended up writing my own digital normalisation method to deal with mapped reads. This method reduces the genome coverage, but preserves reads where coverage is low.

Because I was working with IonTorrent reads (which have variable length), I came up with the idea of selecting the longest read that mapped to each location in the genome (assuming such a read existed). This meant that the highly variable coverage for different samples (sometimes as low as 200X, sometimes as high as 4000X) was flattened out to a much more consistent coverage of about 100-200X. Here's the core of the Perl code that I wrote:

  if(($F[2] ne $seqName) || ($F[3] != $pos) || (length($bestSeq) <= length($F[9]))){
    if(length($bestSeq) == length($F[9])){
      ## reservoir sampling with a reservoir size of 1
      ## See https://en.wikipedia.org/wiki/Reservoir_sampling
      ## * with probability 1/i, keep the new item instead of the current item
      $seenCount++;
      if(!rand($seenCount)){
        ## i.e. if rand($seenCount) == 0, then continue with replacement
        next;
      }
    } else {
      $seenCount = 1;
    }
    if(($F[2] ne $seqName) || ($F[3] != $pos)){
      if($output eq "fastq"){
        printSeq($bestID, $bestSeq, $bestQual);
      } elsif($output eq "sam"){
        print($bestLine);
      }
    }
    $seqName = $F[2];
    $pos = $F[3];
    $bestLine = $line;
    $bestID = $F[0];
    $bestFlags = $F[1];
    $bestSeq = $F[9];
    $bestQual = $F[10];

The full code (which works as a filter on uncompressed SAM files) can be found here.


I wrote a very quick and dirty script to handle conversion between file types using BioJava.

https://github.com/eedlund/Utils/tree/master/BioUtils

Download the jar file here

To run: java -jar BioUtils.jar $FILE $TYPE

where \$FILE is a PDB or mmCIF file you'd like to convert and \$TYPE is the format of the output file [PDB, CIF, MMTF].


Assembly reads having a copy of their beginning in their tail

If you google AGTATGTACAAATACCTACAACTTGTGCT you'll find it's a primer sequence: https://artic.network/resources/ncov/ncov-amplicon-v3.pdf


If nucleotides are acceptable in prosite format because I thought it was protein motifs, the most accurate would be,
A-T-[AT]-G-x-C-[AGC]-C-[GCA](0,1,3)-x-A

The key point is that,

[GCA](0,1,3) of x(0, 1, 3) means zero/ nucleotides (3 indels) or one/ G or C or A nucleotide, and 2 indels or three/ G or C or A nucleotides.
x is any nucleotide i.e. A or G or T or C - it doesn't mean an indel. x(0, 1) is the best representation of a single indel (you call them gaps), here meaning any nucleotide may or may not be present at a given site.

Example from Prosite documentation:

<A-x-[ST](2)-x(0,1)-V
This pattern, which must be in the N-terminal
of the sequence ('<'), is translated as:
Ala-any-[Ser or Thr]-[Ser or Thr]-(any or none)-Val


different gene expression in the monocle

Seurat FindMarker Doc
Yes based upon the R documentation here it seems possible.

This  Differential expression article  seems to show this is how.

library(Seurat)
pbmc <- readRDS(file = "../data/pbmc3k_final.rds")
# list options for groups to perform differential expression on
levels(pbmc)
# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
monocyte.de.markers <- FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono")
# view results
head(monocyte.de.markers)
# Find differentially expressed features between CD14+ Monocytes and all other cells, only
# search for positive markers
monocyte.de.markers <- FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = NULL, only.pos = TRUE)
# view results
head(monocyte.de.markers)
# Pre-filter features that are detected at <50% frequency in either CD14+ Monocytes or FCGR3A+
# Monocytes
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", min.pct = 0.5))
#...
# Test for DE features using the MAST package
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", test.use = "MAST"))
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))

Maybe something like this would work for you.  Monocle export

  #replace the monocle_cds with your monocle
  seurat <-exportCDS(monocle_cds, export_to = c("Seurat", "Scater"))
  #This bellow will list the options for ident.1 and ident.2
  levels(seurat)
  # insert name from levels(seurat) command in parentheses
  head(FindMarkers(seurat, ident.1 = "Double-KO+ Mono", ident.2 = "WT+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))

Hopefully you are able to separate the two groups Double-KO and WT into ident.1 and ident.2 so you can use the function on them. Good luck. As long as the exportCDS function works on your data type you input for cds.


Pathway level analysis of single-cell gene expression

Package Hipathia allows you to compute a value of "activity" for each pathway and each cell, so that the matrix of gene expression is transformed into a matrix of pathway activity. You can use this matrix for further analysis, cell differentiation and clustering. You can also compute functional activity matrices in the same way.


Arbitrary record access in constant time
To get a random record in constant time, it is sufficient to get an arbitrary record in constant time.
I have two solutions here: One with tabix and one with grabix. I think the grabix solution is more elegant, but I am keeping the tabix solution below because tabix is a more mature tool than grabix.
Thanks to user172818 for suggesting grabix.
Update
This answer previously stated that tabix and grabix perform lookups in log(n) time. After taking a closer look at the grabix source code and the tabix paper I am now convinced that lookups are independent of n in complexity. However, both tools use an index that scales in size proportionally to n. So, the loading of the index is order n. However, if we consider the loading of the index as "...a single limited transformation of the data to another file format...", then I think this answer is still a valid one. If more than one record is to be retrieved, then the index needs to be stored in memory, perhaps with a framework such as pysam or htslib.
Using grabix

Compress with bgzip.
Index the file and perform lookups with grabix

In bash:
gzip -dc input.fastq.gz | bgzip -c > output.fastq.gz

grabix index output.fastq.gz

# retrieve 5-th record (1-based) in log(n) time
# requires some math to convert indices (4*4 + 1, 4*4 + 4) = (17, 20)
grabix grab output.fastq.gz 17 20

# Count the number of records for part two of this question
export N_LINES=$(gzip -dc output.fastq.gz | wc -l)

Using tabix
The tabix code is more complicated and relies on the iffy assumption that \t is an acceptable character for replacement of \n in a FASTQ record.  If you are happy with a file format that is close to but not exactly FASTQ, then you could do the following:

Paste each record into a single line.
Add a dummy chromosome and line number as the first and second column.
Compress with bgzip.
Index the file and perform lookups with tabix

Note that we need to remove leading spaces introduced by nl and we need to introduce a dummy chromosome column to keep tabix happy:
gzip -dc input.fastq.gz | paste - - - - | nl | sed 's/^ *//' | sed 's/^/dummy\t/' | bgzip -c > output.fastq.gz
tabix -s 1 -b 2 -e 2 output.fastq.gz

# now retrieve the 5th record (1-based) in log(n) time
tabix output.fastq.gz dummy:5-5

# This command will retrieve the 5th record and convert it record back into FASTQ format
tabix output.fastq.gz dummy:5-5 | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'

# Count the number of records for part two of this question
export N_RECORDS=$(gzip -dc output.fastq.gz | wc -l)

Random record in constant time
Now that we have a way of retrieving an arbitrary record in log(n) time, retrieving a random record is simply a matter of getting a good random number generator and sampling. Here is some example code to do this in python:
Using grabix
# random_read.py
import os
import random

n_records = int(os.environ["N_LINES"]) // 4
rand_record_start = random.randrange(0, n_records) * 4 + 1
rand_record_end = rand_record_start + 3
os.system("grabix grab output.fastq.gz {0} {1}".format(rand_record_start, rand_record_end))

Using tabix
# random_read.py
import os
import random

n_records = int(os.environ["N_RECORDS"])
rand_record_index = random.randrange(0, n_records) + 1
# super ugly, but works...
os.system(
    "tabix output.fastq.gz dummy:{0}-{0} | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'".format(
        rand_record_index)
)

And this works for me:
python3.5 random_read.py

Disclaimer
Please note that os.system calls a system shell and is vulnerable to shell injection vulnerabilities. If you are writing production code, then you probably want to take extra precautions.
Thanks to Chris_Rands for raising this issue.


What is the distribution of indel sizes in a healthy human genome? of insertion:deletion ratios?

One of the 2015 papers from the 1000 genomes project has a nice figure (figure 1) showing the size distribution of medium to large sized insertions and deletions:

From another 2015 1000 genomes paper, one can see that the absolute number of smaller indels is much larger, though an exact size range isn't given (as far as I saw). If you really want to know that, just download the most recent 1000 genomes VCF file(s) and compute the exact median size and/or distribution that you want.