How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?

1. Adapter Trimming

One of the first things I do after encountering a set of reads is to remove the adapter sequences from the start and end of reads. Most basecalling software includes some amount of built-in adapter trimming, but it is almost always the case that some adapter sequence will remain. Removing adapters is helpful for mapping because it reduces the effort the mapper needs to go through to make a match (which may help some borderline reads to be mapped properly).

My preferred adapter trimmer is Trimmomatic. It has the ability to search for palindromic adapter read through, and includes an adaptive sliding window option. Trimmomatic is threaded, and seems to work reasonably fast when run through files one at a time in threaded mode. Here's an example script that I might run to trim reads from some samples. This puts trimmed FASTQ files into a trimmed subdirectory:

TRIMPATH=/data/all/programs/trimmomatic/Trimmomatic-0.36;
JARPATH="${TRIMPATH}/trimmomatic-0.36.jar";
ADAPTfile="${TRIMPATH}/adapters/TruSeq3-PE-2.fa";
CMDS=" ILLUMINACLIP:${ADAPTfile}:2:30:10:5:true LEADING:3 TRAILING:3  SLIDINGWINDOW:10:20 MINLEN:40 ";

mkdir -p trimmed;

for INFILE1 in *_1.fastq.gz; do
  base=$(basename "$INFILE1" _1.fastq.gz);
  echo ${base};
  INFILE2="${base}_2.fastq.gz";
  OUTFILE_P1="trimmed/${base}_P1.fastq.gz";
  OUTFILE_P2="trimmed/${base}_P2.fastq.gz";
  OUTFILE_U1="trimmed/${base}_U1.fastq.gz";
  OUTFILE_U2="trimmed/${base}_U2.fastq.gz";
  java -jar "${JARPATH}" PE -threads 20 -phred33 \
        "${INFILE1}"  "${INFILE2}" \
        "${OUTFILE_P1}" "${OUTFILE_U1}" "${OUTFILE_P2}"  "${OUTFILE_U2}" \
        "${CMDS}";
done;

Trimmomatic is somewhat sensitive to bad input data, and will loudly fail (i.e. stop running) if it sees quality strings of different length to sequence strings:

Exception in thread "Thread-1" java.lang.RuntimeException: Sequence and quality length don't match: 'TACATGGCCCTGAAATGACTTTCACCCAGGCAACCAGTGCCCCCTGTATAGACACATGCCTTGGGCGCTCCCCACCCTTCCTCGCGTGGCCACACCTCTGT' vs '-AAFFJFJJ7A-FJJJFJFJFJJFFA-FFAF-<A-<FFFJA-<-A7-F<<FFJJAJJJJJJJJJ--<--7-7AJ7<AAJA--J7<ACTGCTGTGGGGCACCCAGCCCCCCAGATAGCCTGGCAGAAGGATGGGGGCACAGACTTCCCAGCTGCACGGGAGAGAC'
        at org.usadellab.trimmomatic.fastq.FastqRecord.<init>(FastqRecord.java:25)
        at org.usadellab.trimmomatic.fastq.FastqParser.parseOne(FastqParser.java:89)
        at org.usadellab.trimmomatic.fastq.FastqParser.next(FastqParser.java:179)
        at org.usadellab.trimmomatic.threading.ParserWorker.run(ParserWorker.java:42)
        at java.lang.Thread.run(Thread.java:745)

2. Read Mapping

Read mapping finds the most likely location for a read within a target genome.
There are extremely fast approximate mappers available, but these don't yet work for variant calling, and an exact mapping approach is necessary. My current preferred mapper is HISAT2. I use this instead of BWA due to the double-read issue, and because of the local variant-aware mapping. HISAT2 is made by the same computing group as Bowtie2 and Tophat2 (JHUCCB), and is their recommended tool for replacing those other programs.

  HISAT2 is a successor to both HISAT and TopHat2. We recommend that HISAT and TopHat2 users switch to HISAT2.

The HISAT2 page includes a link to a genomic index file for the human genome, including SNP variants, which I use when mapping human reads. HISAT2 is also threaded, so I run it on the files one at a time and pipe through samtools to create a sorted BAM file:

mkdir mapped
for r1 in trimmed/*_P1.fastq.gz;
  do base=$(basename "${x}" _P1.fastq.gz);
  r2="trimmed/${base}_P2.fastq.gz";
  echo ${base};
  hisat2 -p 20 -t -x /data/all/genomes/hsap/hisat_grch38_snp/genome_snp -1 \
    "${r1}" -2 "${r2}" 2>"mapped/hisat2_${y}_vs_grch38_stderr.txt" | \
    samtools sort > "mapped/hisat2_${y}_vs_grch38.bam";
done

3. Variant Calling

I currently use samtools/bcftools to do this, but would be interested in other opinions as we recently had a grant reviewer response that samtools was a dinosaur program and better approaches were available. Samtools doesn't currently work with threads for variant calling, so I save the commands to a text file and then run them through GNU Parallel. This requires the genome files to be first downloaded and indexed via samtools faidx. This produces a set of gzipped VCF files in the variants directory:

mkdir -p variants;
(for x in mapped/*.bam;
   do echo samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa "${x}" \| \
   bcftools call --ploidy GRCh38 -v -m -O z -o variants/$(basename "${x}" .bam).vcf.gz; done) > call_jobs.txt
cat call_jobs.txt | parallel -u;

As a slow, but more accurate alternative, variants can be called for all samples at the same time:

samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa mapped/*.bam | \
  bcftools call --ploidy GRCh38 -v -m -O z -o variants/hisat2_allCalled_vs_grch38.vcf.gz

If a faster processing is desired, then the variant calling can be parallelised across chromosomes using parallel and merged afterwards using bcftools norm. The implementation of this is left as an exercise to the reader.


Why is library(Seurat) causing my R session to have a fatal error?

If you haven't already done so I would make a complete installation of Anaconda, particularly good with bash. The I would either install using Anaconda Navigator, or ...

# Create an R environment and activate it
conda create -n r_env r-essentials r-base python=3.7
conda activate r_env

# install Reurat
conda install -c bioconda r-seurat

# installing standard packages
conda install -c r r-randomforest --yes

If you have a Mac the instructions are here. In some systems you need to just activate r_env .

Please note if you've installed e.g. python3 with pip a conflict occurs with the Anaconda package.


The problem is the scale used: For the plot you called "weird" (first from the top), the scale is 50 and for the "ggplot only" (third from the top) the scale is 1.

You should play with the stat_compare_means(label.y = 50) bit, you can try setting the label.y parameter to 1.5 or 2.

EDIT BECAUSE OF THE ADDITIONAL REMARKS ON THE QUESTION:

Once again (and for sure the last time), I had to manually create example data as OP ignores repeated reminders on how to make a reproducible example with R.

library(data.table)
library(ggpubr)

text <- "
treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
1       post     Sig1 0.1173863
2       post     Sig1 0.1888243
3       post     Sig1 0.2191765
4       post     Sig1 0.1277735
5       post     Sig1 0.1711712
6       post     Sig1 0.0000000
1021      pre     Sig9     0
1022      pre     Sig9     0
1023      pre     Sig9     0
1024      pre     Sig9     0
1025      pre     Sig9     0
1026      pre     Sig9     0
1021      post     Sig9     1
1022      post     Sig9     1
1023      post     Sig9     1
1024      post     Sig9     1
1025      post     Sig9     1
1026      post     Sig9     1
"
my_data <-fread(text)

> my_data
      V1 treatment variable     value
 1:    1       pre     Sig1 0.1173863
 2:    2       pre     Sig1 0.1888243
 3:    3       pre     Sig1 0.2191765
 4:    4       pre     Sig1 0.1277735
 5:    5       pre     Sig1 0.1711712
 6:    6       pre     Sig1 0.0000000
 7:    1      post     Sig1 0.1173863
 8:    2      post     Sig1 0.1888243
 9:    3      post     Sig1 0.2191765
10:    4      post     Sig1 0.1277735
11:    5      post     Sig1 0.1711712
12:    6      post     Sig1 0.0000000
13: 1021       pre     Sig9 0.0000000
14: 1022       pre     Sig9 0.0000000
15: 1023       pre     Sig9 0.0000000
16: 1024       pre     Sig9 0.0000000
17: 1025       pre     Sig9 0.0000000
18: 1026       pre     Sig9 0.0000000
19: 1021      post     Sig9 1.0000000
20: 1022      post     Sig9 1.0000000
21: 1023      post     Sig9 1.0000000
22: 1024      post     Sig9 1.0000000
23: 1025      post     Sig9 1.0000000
24: 1026      post     Sig9 1.0000000
      V1 treatment variable     value

ggboxplot(my_data,
          x = "treatment",
          y = "value",
          color = "treatment",
          facet.by = "variable") +
  stat_compare_means(label = "p.format")

And this code generates:


For phylogenetic tree construction from core-genome which one is preferable: amino-acid based MSA or nucleotide based MSA?

For a species complex nucleotide phylogenies all the way.

The reason is neutral mutation, which is not observed at the amino acid level because these mutations are part of the degenerate code which don't induce amino acid substitutions. Such mutations are commonly referred to as the third codon wobble, although occasionally the first codon can mutate and not produce a change in amino acid.

Neutral otherwise known as synonymous mutations occur at a much higher rate than mutations that cause a substitution of an amino acid. This is the reason why  you were told information is lost at the amino acid level because you're losing neutral mutations.

There are circumstances when the amino acid phylogenies are preferred to nucleoside trees, but these are bit pretty specialised areas.


Retrieving NCBI Taxa IDs from refseq or GenBank assembly accession

Turns out I'd already written some code that did it, human memory is a funny thing.

This uses biopython to split the field description to where the species is. May not work for all NCBI files, but seems to work on most.

import Bio
from Bio import SeqIO
from Bio import AlignIO

for record in SeqIO.parse (FILE, "fasta"):
    Speciesname = record.description.split('[', 1)[1].split(']', 1)[0]


How can I install bioconductor-gviz and use it in jupyter notebook?

# Install the BiocManager package first
install.packages("BiocManager")
library("BiocManager")

# Then install the Gviz package
BiocManager::install("Gviz")


Kruskal-Wallis rank sum test output confusion

Your problem comes from using an incorrect test, the formula parameter says:

  a formula of the form response ~ group where response gives the data values and group a vector or factor of the corresponding groups.

When you apply the kruskal  test, you don't have groups in the LSC variables! You aren't comparing the mean of several groups based on LSC!
When you use stat_compare_means it is doing a wilcox.test (it hints to it in the help page "a list of additional arguments used for the test method. For example one might use method.args = list(alternative = "greater") for wilcoxon test.").

When you use the wilcox.test we get the same p-value:

wilcox.test(df1$HSC, df1$LSC)

    Wilcoxon rank sum test with continuity correction

data:  df1$HSC and df1$LSC
W = 110980, p-value = 4.849e-12
alternative hypothesis: true location shift is not equal to 0

In addition, probably you should read about the paired argument, if this two values where measured on the same samples you should use it (but I can't tell it for sure because I don't know what does it mean each variable name).


Drawbacks of upper quartile normalization for scRNA-seq data

My primary concern with using the top ~1% or so in upper quantile normalization is that it's going to be prone to the same robustness issues that RPKM/FPKMs run in to. That is, if for whatever technical reason you have to have a fair bit of variability in a couple very highly expressed genes (typically rRNAs, but one can imagine other genes) and the set of genes going into the UQ normalization is composed primarily of those then the normalization results will simply reflect irrelevant technical variability.

Of course, this all assumes that it's even appropriate to normalize scRNAseq data with a standard method like this. In some cases it is (e.g., when comparing control vs. treatment with a largely homogeneous cell population), but if you're dealing with very different cell types then I worry that the standard method is going to wash away things like differences in bulk RNA content between cell types. Presuming there are spike-ins, they should have the same normalization applied and then checked to ensure there's nothing greatly amiss. Having said that, if you're going to be doing tSNE and want to look for cell types, then odds are good you'll still be able to see most of them (after all, it's unlikely that RNA quantity is the only difference between them, though not removing that difference would make things easier).

The problem of global shifts in expression

The concept of global shifts in expression is not unique to scRNAseq, in fact the most common example of this is Myc up/down regulation in relation to cancer. To explain the issue, let's suppose that you have the following (very small) matrix of counts:

       sample1 sample2
gene1        1       2
gene2       10      20
gene3       15      32
gene4      100     198

As you can see, sample2 has roughly twice as many reads for each gene as sample1, so any normalization is going to cut all of its counts in half (or the reverse for sample1, or perhaps a mixture of the two...depending on the method). That's usually fine, since the assumption with all common normalization methods is that differences like this are simply due to differences in sequencing depth between the samples. But what happens if sample2 actually has twice as much RNA as sample1? Well, then you've just normalized away all of the differences between the samples (I encourage readers to go through the literature on Myc, since that's a very nice example of this, where a mostly-global shift in expression led to apparent down-regulation of genes with constant expression). Such a situation doesn't occur very commonly in bulk RNAseq, but when it does you need to use something like spike-ins (dosed to match cell number or DNA quantity) to maintain the global difference.

If this is only rarely an issue with bulk RNAseq, then why is it more likely to be an issue in scRNAseq? The answer is that bulk RNAseq benefits from the averaging out of many different cell-types, such that ensemble averages/medians don't tend to change much due to treatments...but this isn't going to be the case at the single-cell level, since cells can vary vastly in size (compare a sperm and an oocyte...every RNA should be differentially expressed between them). This difference in base RNA content between different cell types will be washed away by common normalization methods, which can lead to either losing differences or creating differences in genes where there shouldn't be (again, see examples of this in the Myc literature). Whether this really becomes an issue will depend on what your goals are with your data and how much of a global expression difference you actually have and the type of analysis you want to perform.

At the end of the day, just remember that common RNAseq normalization methods assume that (A) there's no global expression change between samples and that (B) if there are very many DE genes that the direction of change is not heavily skewed in one direction. If either of these assumptions are violated then you can run into issues.


How to append numbers only on duplicates sequence names?

And here’s a solution using R (with the Bioconductor):

fa = ShortRead::readFasta(infile)
ids = as.character(ShortRead::id(fa))
fa@id = Biostrings::BStringSet(make.unique(ids, sep = '_'))
ShortRead::writeFasta(fa, outfile)


I think there needs to be clarification of the question. Shannon entropy is not the same thing as the Shannon-Wiener diversity index (H), commonly referred to as the Shannon index.

The reason I brought it up is because you have the tag "population genetics". Shannon diversity index is commonly used in population genetics to compare species richness (the number of species in a given area) and their relative abundance.

Shannon entropy is more to do protein structure, which isn't really population genetics, but can relate to how stable a given mutation might be. That is not really population genetics however, population genetics is about abundance of a given allele in this context.

I am pretty sure that on this site most responders will recommend alternative free energy approaches to perform this style of analysis.


the fact is that the present upgma tree as a midpoint rooting here, corresponding to node 46.

getRoot(upgma.a)

so for this kind of tree, I recommend you to use the midpoint() rooting function.

pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
plot(midpoint(pars.a))


Hba-a2 has two copies on UCSC genome browser

Gencode is correct (when in doubt, always assume Gencode/Ensembl is correct when it differs from UCSC). What presumably happened is that the UCSC annotation pipeline noticed that Hba-a1 and Hba-a2 are nearly identical (they produce the same protein and have minimal nucleotide difference) and just misannotated.


what does a sequence look like before alignment?

Sequence alignment typically refers to the primary sequence. That is, we ignore the complex 3D structure of the molecule and its spatial arrangement in the cell. Instead, we use a symbolic representation composed of symbols like A, C, G, T, and U. Alignment is a conceptual tool that, as you pointed out, is used to find regions of similarity between sequences. Therefore, it only makes sense in the context of two or more similar sequences.

By way of example, consider a simple case of two sequences X and Y. This is what they would look like before alignment. The similarities between the two sequences are unclear.

X: ATTGACATAGGCACTTCGACTA
Y: GACATCTTCGGCTAA

This is what they would look like after alignment. The similarities (and differences) are now much easier to see.

X: ATTGACATAGGCACTTCCACTA-
      |||||     ||||  |||
Y: ---GACAT-----CTTCGGCTAA

This arrangement identifies regions of similarity that may be related through evolution. This doesn't always mean perfect agreement though. As this example shows, sometimes mismatches or gaps are allowed so that long stretches of identical matches are lined up correctly. Sometimes, for very closely related sequences, the optimal alignment doesn't require any gaps at all. This seems to be the case with some of your recent questions about sequence logos.

Note: Sequence alignment is a fundamental concept in bioinformatics, and there is a tremendous amount of research and development on the topic. It's hard to imagine doing any amount of serious DNA sequence analysis responsibly without a basic understanding of sequence alignment. I would recommend searching Google and/or Google Scholar for intro-level reviews of the topic.


Merging several fq.gz files or R1 and R2 classes into a single one

What you did should work.  You don't need to decompress and recompress.
I'd check the files.  Were they trimmed?  Can you get the top, say 1000 lines of the merged files to work?  (You will need to decompress to count lines accurately)


Nextflow error: can not invoke method minus () on null object

Channel names must consist only alphanumeric characters including underscores. Your channel called 'plink-data' contains a hyphen. This is problematic because it is interpreted as a minus() operation to remove 'data' from 'plink'. You get a NullPointerException since 'plink' hasn't been defined. You can test this easily with the following example:
Channel
    .of(1..23, 'X', 'Y')
    .set { chr-names }

chr-names.view()

Results:
Cannot invoke method minus() on null object

The solution is to just use an underscore instead of a hyphen:
Channel
    .of(1..23, 'X', 'Y')
    .set { chr_names }

chr_names.view()

Follow-up issues:
I think if you want to use the bfiles plink option (and I think this is reasonable), I think it makes the most sense to just keep your bfiles in a tuple:
Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size: 3 )
    .set { plink_data }

process test1{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bfiles) from plink_data

    output:
    tuple val("${sample_name}-t1"), path("${sample_name}-t1.{bed,bim,fam}") into test1_results

    """
    plink \\
        --bfile "${sample_name}" \\
        --maf 0.05 \\
        --make-bed \\
        --out "${sample_name}-t1"
    """
}

process test2{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bfiles) from test1_results

    output:
    tuple val("${sample_name}-t2"), path("${sample_name}-t2.{bed,bim,fam}") into test2_results

    """
    plink \\
        --bfile "${sample_name}" \\
        --maf 0.01 \\
        --make-bed \\
        --out "${sample_name}-t2"
    """
}

If you need to reference each bfile separately, you could just unpack the tuple and reference each file accordingly:
Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size: 3 )
    .map { group_key, files -> tuple( group_key, *files ) }
    .set { plink_data }

process test1{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bed), path(bim), path(fam) from plink_data

    output:
    tuple val("${sample_name}-t1"), path("${sample_name}-t1.{bed,bim,fam}") into test1_results

    """
    plink \\
        --bed "${bed}" \\
        --bim "${bim}" \\
        --fam "${fam}" \\
        --maf 0.05 \\
        --make-bed \\
        --out "${sample_name}-t1"
    """
}

process test2{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bed), path(bim), path(fam) from test1_results

    output:
    tuple val("${sample_name}-t2"), path("${sample_name}-t2.{bed,bim,fam}") into test2_results

    """
    plink \\
        --bed "${bed}" \\
        --bim "${bim}" \\
        --fam "${fam}" \\
        --maf 0.01 \\
        --make-bed \\
        --out "${sample_name}-t2"
    """
}

Note that in the first example, bfiles could be replaced with "*" since no reference is made to the tuple of bfiles anyway. Also note that the list of files produced using the fromFilesPairs operator are returned in lexicographic order.


Splitting into multiple files and changing the IDs can be easily done:

perl -pe 'if(/>/){/\[(.*?)\]\s*$/; $_="> $1\n"}' file.fa |
    awk '(/^>/){name=$2} {print >> name".fa"}'

That assumes all your FASTA headers have [foo bar baz] as the last element of a line. It will create a file called foo.fa (the bacterium's name) with all sequences saved there.

Explanation

perl -pe : run the script given by -e on each line of the input file, and print the resulting line.
if(/>/) : if this line starts with a >. \
/\[(.*?)\]\s*$/ : match an opening bracket (\[), then capture (that's what the parentheses do, they capture a pattern so we can refer to it as $1) everything until the first ] (.*?\])
$_="> $1\n" : the $_ special variable in Perl is (in this case) the current line. So, $_=foo means "make the current line read foo. Since the -p prints each input line, changing the value of $_ means the changed value will be printed. So here, we are printing >, whatever was in the square brackets ($1) and a newline character.

The output of the perl command alone on your example input file is:

$ perl -pe 'if(/>/){/\[(.*)\]\s*$/; $_="> $1\n"}' file.fa
> Arthrobacter phage Mudcat
MGLSNTATPLYYGQF...

> Achromobacter phage phiAxp-3
MSNVLLKQELDEWL...

> Delftia phage RG-2014
MSEPRKLVKKTLD...

So, we now pass it through awk which does:

(/^>/){a=$2} : if this line starts with an >, save the second field (the bacterial species) as the variable name.
{print >> name".fa"} : print each line into a file whose name is the current value of the variable name with a .fa. extension.

If you prefer python scripts to the one-liner approach, you can do the same thing with:

#!/usr/bin/env python
import re
outFile = None
with open("file.fa", "r") as inFile, open(species + ".fa", 'a') as outFile:
    for line in inFile:
        line = line.rstrip()
        if line.startswith('>'):
            regex = re.compile('.*\[((.+?)\s+.*?)\].*')
            matches = regex.search(line)
            species = matches[2]
            outFile.write('>%s\n' % matches[1])
        else:
            outFile.write("%s\n" % line)

As for your script, you've got the right idea, but have a small bug. You forgot to remove the \n from your input file, so it looks for Arthrobacter\n instead of Arthrobacter. The golden rule of debugging is "print all the things". If you add print("Searchterm: ",searchterm)  you will see:

Searchterm:  (terminase large subunit AND viruses[Organism]) ANDArthrobacter
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDAchromobacter
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDDelftia
AND refseq[Filter]

So, you need to remove the newline characters and add a space like so (I also made it a bit more "pythonic" and conforming to the Python syntax guidelines):

#!/usr/bin/env python
from Bio import Entrez
Entrez.email = "example@example.org"

with open("species2.txt", "r") as in_handle:
    for line in in_handle:
        line = line.rstrip()
        searchterm = ("(terminase large subunit AND viruses[Organism]) " +
                      "AND %s AND refseq[Filter]" % line)
        print("Searchterm: ", searchterm)
        searchResultHandle = Entrez.esearch(db="protein",
                                            term=searchterm, retmax=1000)
        searchResult = Entrez.read(searchResultHandle)
        ids = searchResult["IdList"]

        handle = Entrez.efetch(db="protein", id=ids,
                               rettype="fasta", retmode="text")
        record = handle.read()

        with open('terminase_large_' + str(line[:-1]) + '.fasta', 'w') as out_handle:
            out_handle.write(record.rstrip('\n'))


Help interpreting my phylogenetic tree construction of bacterial species in same genus

The top phylogeny is likely mid-point rooted, but I suspect the rooting is fine. The lower diagram has been rooted on the query sequence. It is taxonomic analysis using 16S which in bacteria is considered to represent species delineation at a given percentage of divergence. I thought it was 5%, which is impossible based on the phylogram in the top picture.

Bootstrap values lower than 80% mean that the clade is not supported (monophyletic group is a better term). The delineation between 10001 and 10013 against the massive ingroup of 13 16S sequences is robustly supported by bootstraps of 95% and 99%. The divergence between them is around 2%. There is no taxonomic distinction between the 13 16S sequences based on the phylogram.

One question is whether this is full length 16S? If it is there appears to be a species designation confusion in the lower clade supporting 10001 and 10013 because it is identical to to Dermococcus species and there is no way to assess which one is valid, nor whether they are really members of the same species. The very small divergence suggesting in D. abyssis might also call this into question. 10013 is in fact a robust outgroup to 10001 and the three Dermococcus spp. but the genetic distance involved is small.

However, this is not the major result. You have recovered what appears to be a large population (if that is how you sampled it) with close genetic homology to Dermococcus strain DSM20448. Given the complexity of species lineation within this genus it is impossible to say whether this represents the same species as DSM20448 but you could call it a "species complex".

Overall the population structure (if that is how it was sampled) is heavily skewed towards a bacterial species which is hitherto received minimal attention. This population skew (if that is how it was sampled), could be reversed in different ecotopes/niches/ecoregions/patients?

In summary, I wouldn't get caught up in systematics because in this genus its a bit messy. I would focus on the population structure and these results are very clear.


NaN values after ComBat analysis on TCGA COAD RNAseq

I think I might know what the problem is.

I plotted the log2 corrected values, and the distribution may not be normal because of the 0 values present in RNASeq data. If I understand correctly, this data is actually missing, whereas with my current analysis it is viewed by R as an actual value of 0. So either the missing data needs to be imputed, or the data needs to be scaled and normalized to mean of 0, then subsequently log2 transformed.

> data <- standardNormalization(data)
> data <- log2(data + 1)
> combat_data = ComBat(dat=data, batch=batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE)
> combat_data[1:5,1:2]
          X01240896.3f3f.4bf9.9799.55c87bfacf36
100507661                           -0.13857428
57103                               -0.06700506
22838                                0.16295096
55567                               -0.13579991
6147                                 1.64990772
          X01ad5016.f691.4bca.82a0.910429d8d25b
100507661                           -0.13789534
57103                               -0.06210232
22838                                0.01546224
55567                               -0.13443289
6147                                 1.48473072


vcf-merge (or better: bcftools merge) does a column-based merge of VCF files, basically taking on an additional genotype column to each variant. While this sort-of works, it doesn't take into account all the information about variation. There might be low coverage at a region in one VCF file that indicates a deletion, but the deletion is only observable when considering other samples with high coverage at that same region. It's also possible that INDELs could be defined differently in different samples, requiring an additional bcftools norm step for left-normalisation prior to the merge.

Assuming this column-based merge is what you want (and if not, then bcftools concat is your friend for row-based merging), then it's better to regenerate the merged VCF file as a multi-column file directly. One way of doing this is feeding all BAM files into samtools mpileup to generate genotype likelihoods, then feeding the result of that into bcftools call. Something like the following:

samtools mpileup -uf ref.fa in1.bam in2.bam in3.bam | \
  bcftools call -mv -O z -o variants_ref_1-3.vcf.gz


How to exclude the repetition of gene-gene correlation calculation in python?

You don't need g2 to go from 0 to din.index, just from 0 to g1 - 1. This way, you'll end up calculating just for the "lower triangle":

g1-> 0   1   2   3
g2
0                    #(0,0 to 0,1-1 = nothing)
1    X               #(1,0 to 1,1-1 = 1,0)
2    X   X           #(2,0 to 2,2-1 = 2,0; 2,1)
3    X   X   X       #(3,0 to 3,3-1 = 3,0; 3,1 and 3,2)
4    X   X   X   X   #(4,0 to 4,4-1 = 4,0; 4,1; 4,2 and 4,3)


Why do I have >10,000 cells in the 10X matrix produced by cellranger?

The raw data from cell ranger contains all of the barcodes detected in the experiment.  These raw barcodes are filtered by cell ranger to identify the barcodes that likely represent cells rather than empty droplets / dead cells. You can find the filtered data in the filtered feature-barcode matrices.  You can also examine the Run summary HTML file to see a summary of the cell ranger barcode filtering.

The reason why there are so many raw barcodes is an important optimization by 10X to increase the cell capture rate compared to techniques such as drop-seq.  In a 10X experiment the number of gel beads is much higher than the number of cells which leads to a higher probability of a cell and a bead ending up in the same droplet.  This also leads to a large number of empty droplets with just a bead and some ambient RNA or a dead cell.  These all end up being sequenced and the empty drops tend to have a small number of unique molecules (the majority less than 10 UMIs / barcode).


Problem installing last build of Sniffles

You need to install GCC, or a newer version of Clang.

Your compiler — though the binary is called gcc — is actually Clang, and the Clang version that’s shipping with macOS/Xcode still does not support OpenMP.

The easiest course of action is probably to install GCC via Homebrew (brew install gcc) and then configure Sniffles to use Homebrew’s GCC (i.e. setting the CMAKE_C*_COMPILER paths to /usr/local/bin/g{cc,++}-$version).


I'm pretty sure there are no public databases of variants that aren't included in dbSNP. That is the variant database so that's the one you want to use.

While there are various other databases around, they tend to use dbSNP rsIDs as well so any variants in those will also be in dbSNP. So while dbSNP is nowhere near complete, of course, it's still the most complete one you're likely to find.


There is a conda package for maxquant that should work on most platforms (I've used it on CentOS and I think it's been used on Ubuntu as well). There's also a docker container version of that that we (Bioconda) have posted to the biocontainers project that should also work.

Please note that large parts of a maxquant analysis are single threaded, regardless of what you specify. We've also had single analysis runs lasting over a month (even when using 90 threads), so the program can be very slow if your search space is large.


How are centromeric/telomeric regions categorized in RepeatMasker? "Other" category?

On (2), centromeres and telomeres primarily consist of satellite sequences. In human centromeres, the dominant type is alpha satellite (aka alphoid). GGAAT micro-satellite (aka SATIII centromeric micro-satellite) is also very common. There are less common types such as beta satellite. Because there are all kinds of centromeric/telomeric sequences of different lengths and properties, RepeatMasker doesn't put them into one class. Note that RepeatMasker has a "centr" class but it does not include all centromeric satellites (no SATIII, for example). In general, it is not advised to find centromeres with RepeatMasker output.

GRCh38 is an unusual reference genome in that its centromeric sequences are computationally generated modeled after Venter as I remember. They are mostly alpha satellite, but are not real sequences (which has always bugged me). GRC know exactly where they put these sequences. GRCh37 is a more typical reference genome. It lacks probably majority of alpha satellites, but still has many. When you map reads to GRCh37, you occasionally see ~10,000-fold read depth even though the average depth is ~30X. These regions are mostly alpha satellites. Mappers just pile all different copies of satellites onto the ones remained in GRCh37.


Where can I find a database that has phenotype information together with associated SNPs?

International Mouse Phenotyping Consortium is building a database of phenotypes and knock-outs of mouse. I believe that this database will be fairly complete (20000 knock-outs), but these are knock-outs, not SNPs... There are several mouse GWAS studies, but I am not aware of a database that would pull all the results together.

Arabidopsis big GWAS project is called 1001 genomes project. A lot more information about Arabidopsis from individual experiments is collected and maintained by The Arabidopsis Information Resource, however this resource is subscription based, therefore what you need might not be available for free.

As mentioned, SNPs of rice were identified and linked to different rice groups in The 3,000 rice genomes project. There is also GWAS done on 512 individuals, not sure if it is sufficient for reliable conclusions (even they claim that they found supersigificant links).

Is this what you were searching for?


Mapping RefSeq-ID or NCBI-ID (Entrez-ID) to EC number

You can use Entrez Direct for this as follows:

efetch -db protein -id YP_805528.1 -format gpc -mode xml | xtract -insd Protein EC_number
YP_805528.1     6.3.2.13


How can I predict if the binding site of a DNA-binding protein is located on the promoter sequence of a gene?

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

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

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

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

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

You can define your promoter regions however you like, e.g., from the TSS to 2.5kb upstream of the TSS. It's up to you and your experiment.

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


Fix dash/empty ALT alleles for deletions in VCF

According to the specs the ALT is not allowed to contain -.
Here is a little python (>=3.6) script that is able to fix it. It makes use of the pysam module. You need a reference file and this fasta file must be indexed by samtools faidx reference.fa. Please note that it doesn't support multiple values in the ALT column for records that contain a -.
Usage:
python fixvcf.py -f reference.fasta -o out.vcf in.vcf

The script:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import sys

import pysam

class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """
        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) -> pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """
        if "-" in record.alts:
            ref = self.fasta.fetch(region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1})")
            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "-" else f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see: https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """
        print(self.vcf.header, end="")
        print(*map(self.fix_record, self.vcf), sep="", end="")

if __name__ == "__main__":
    app = Application()
    app.start()

A lot of the code above is boilerplate for running a script on the command line. The critical code for the transformation is in the fix_record function. The relevant code is below, expanded, with line numbers and a description of each command.
[1]   ref = self.fasta.fetch(
          region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1})"
      )
[2]   record.pos = record.pos - 1
[3]   record.ref = f"{ref}{record.ref}"
[4]   record.alts = tuple(
          map(
              lambda alt: ref if alt == "-" else f"{ref}{alt}",
              record.alts
          )
      )

This grabs the base preceding the inserted/deleted sequence.
This corrects the position of the variant to account for prepending another nucleotide.
This defines the new corrected REF allele: concatenate the base grabbed in [1] with the original REF allele
Apply an anonymous (lambda) function to each ALT allele; the function will replace - with the base grabbed in [1], or preprend the base to the allele if it is not -


Problems in creating desired phylogenetic tree with ggtree


How to index residues by icode in biopython?

Okay, I found the answer myself. It is possible to index not only by RES_NUM, but by its id (check with residue.get_id()), so the below example would work:
#structure[0][RES_CHAIN][het, resseq, icode]
structure[0][RES_CHAIN][' ', 27, 'A']


How to get rows with similar values in two different columns using command line?

I think it is easy to do with awk, e.g. see here.
For your example, with your data in file:
%  awk '$2 != $3 {print $0}' file
# outputs
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

%  awk '$2 == $3 {print $0}' file
# outputs:
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Note that where you say "similar" I am interpreting that to mean "exactly the same", because that is the case in your example.


Things like this might depend on your specific library prep, but in general: sequencing starts at the end of the adapter, not before it. You will only see adapters if you sequence through the entire fragment into the adapters on the other side. If your fragment is long enough you won't see any adapter, and of course, that's the desirable outcome of your library prep.


mapping and de novo assembly of plant without reference genome

You have multiple questions, I will try to answer one by one:

On creating an index with STAR: As long as you provide a fasta file and its corresponding gff, you will be able to generate an index.
On creating a reference with Trinity: You have mentioned "to construct the Reference Genome with Trinity", however, with your RNA-seq data and Trinity, you can construct a de novo transcriptome and not a genome.
Uploading files to Galaxy: I suggest to do so via FTP using FileZilla or something similar. The transfer can take some time if your files are large though.

If I understand your questions correctly, you are performing RNA-seq and you are concerned about the completeness of the genome of your model organism. In that case, I would advise to bring together a de novo transcriptome from your RNA-seq data (with Trinity), use one of the quasi-mapping methods (with a tool like Salmon or kallisto).


There are actually two questions.

The first one is about the memory issue. I believe merge from data.table would solve that issue.

The second one would be aggregating or summarizing the identifiers. As @StupidWolf mentioned, you need to have a rule, summing them up, averaging them, ... For this second part one way is using dplyr:

merged_data %>%
  group_by(probe_colum) %>%
  summarise_all(mean or sum or ...) %>%
  mutate(new_column = paste(probe_column, gene_column, sep = "_"))

The same thing can be done faster with data.table:

merged_data[, lapply(.SD, mean or sum or ..., na.rm=TRUE), by = probe_column ] [, new_column := paste(probe_column, gene_column, sep = "_"),]

.SD means subset of data and your subsets are your groups based on your probes.


Error in as.vector(x) : no method for coercing this S4 class to a vector

You can't cbind a bunch of obscure object types.

If you want merged count tables you should do this:

mdat <- do.call(cbind,lapply(dat,assay))

Where row.names are Ensembl gene IDs and col.names are the SRR accessions.

Then run your table writing command.

-

If you want the coordinates of your genes then do this to make a bed with the genomic locations:

library(rtracklayer)
export.bed(rowRanges(dat[[1]]),'gene_coords.bed')


De novo motif discovery in protein sequences

I think MEME is a good tool for your purpose, but there are others as well. I can think of InterProScan for example (although I am not sure if it is really de novo). Here a summary of available tools.


How to deal with correlated SNPs in GWAS?

What you are attempting to do is known as LD-pruning.
As @Emily_Ensembl said, it is not customary to do this for standard association tests: it is possible that one of the SNPs you remove is causal, or a better proxy for the causal locus, and would give (slightly) better association signal than the other. Even for SNPs in perfect LD, pruning is unwise because it would complicate interpretation: if SNP chr1_2 had some obvious function, but you removed it and are left only with its proxy chr1_1 (which could be far away), you will have a hard time connecting the dots. Typical GWAS tests are handled pretty quickly on modern computers, and are linear in SNPs, so there is no computational reason to prune either.

On the other hand, many fancier procedures would require LD-pruned data, because they either
a) are very computationally expensive, or
b) assume independent SNPs.

I suggest reading PLINK's LD pruning documentation, as it pretty much covers the standard ways to do this.


How can I find the position of every mutation where the Allele Frequency is greater than X, in regards to a reference such as Hg19?

Hello and welcome to the world of bioinformatics,

BAM files are compressed formats of SAM and you can read more about it here.

To answer your question, I presume that your BAM is generated from a Whole Genome library preparation. Calling mutations on a Whole Genome takes a longer time than a Whole Exome since you are looking at every position instead of the exons of your genome.

What you are looking for is a special format of files called Variant Caller Format (VCF) generated by Variant Callers tools listed below in comparison studies. There isn't a tool better than another, they all have their advantages and disadvantages. I suggest that you read about them as they also describe their pipeline approaches.

Comparison of three variant callers for human whole genome sequencing
Comparing the performance of selected variant callers using synthetic data and genome segmentation
Performance assessment of variant calling pipelines using human whole exome sequencing and simulated data
Systematic comparison of germline variant calling pipelines cross multiple next-generation sequencers
Variant Callers for Next-Generation Sequencing Data: A Comparison Study

Once the VCF is generated, you will also have a column named AF which corresponds to the allele frequency also computed when using these tools. There are tools that allow you to filter VCF files to extract the needed information with a wished criteria like SnpSift,  vcftools, VariantFiltration. So, I suggest reading the above articles as a start.

As for building pipelines, it depends on what language you feel at ease at. However, if I could suggest, Nextflow is a good start. You can also use Snakemake if you like Python. Or maybe simply bash ? In any case, there are many pipelines already made to help you through your needs.

(If this helps, don't forget to either accept the answer, vote the answer or both).


Get results of keyword search on Pfam via python script

If I understand your question correctly, then you can use the Pfam proteom search to get the data you are looking for.

Go to the

  Pfam website -> Browse -> select the letter D from proteoms category
  -> D. rerio

There you can find a download link under the section Domain composition:

  You can download a list of all regions for this proteome from our FTP
  site.

The download link seems to be bound to the species ID and if you want to automate your data retrieval you may be able to simply download the file from there in your python script.


How can I calculate gene_length for RPKM calculation from counts data?

Here you can find some example R code to compute the gene length given a GTF file (it computes GC content too, which you don't need). This uses one of a number of ways of computing gene length, in this case the length of the "union gene model". In this method, the non-duplicated exons for each gene are simply summed up ("non-duplicated" in that no genomic base is double counted). This is a very simple way of getting a gene length.

There are alternative methods that you should be aware of, among which are:

Median transcript length: That is, the exonic lengths in each transcript are summed and the median across transcripts is used. This is probably a little more valid than the code that I linked to.
Per-sample effective gene lengths: the optimal method, though it requires using something like RSEM, which will give you an effective gene length.
Using the length of the "major isoform" in your tissue of interest. This isn't as good as method 2, but is more accurate than all of the others.

At the end of the day, you're just coming up with a scale factor for each gene, so unless you intend to compare values across genes (this is problematic to begin with) then it's questionable if using some of the more correct but also more time-involved methods are really getting you anything.

Edit: Note that if you want to plug these values into some sort of subtyping tool (TNBC in your case), you should first start with some samples for which you know the subtype. Then you can at least see if you're getting reasonable results. After that, do read up on how the method works and see if there's anything about RNAseq that makes it incompatible.


I guess you already know what identity means. But for completeness, it is the % of pairwise alignment position which have an identical amino acid residue. For the identity+similarity it is important to know that amino acids can be grouped according to their physicochemical properties e.g. charge, size etc. While there are standard groupings for amino acids it may be that the proprietary algorithm uses different groupings. Either way, % identity+similarity would imply the number of pairwise alignment positions which contain either the same amino acid or amino acids from the same group (however they define groups).


I'm not a fan personally, but it is context specific. The problem is that it is heading towards optimisation, particularly algorithmic optimisation. If the alignment is the final calculation - that is cool, but the situations are select. It depends how far you want to take it.
The rationale of my skepticism is simple - if an alignment was decided on a 'maximization' criteria, be it likelihood or some other metric - the subsequent calculations based on the alignment, such as phylogenetics Jukes Cantor, Kimura 2 and 3-parameter etc ... we enter circular mathematics. The critiera that defined the alignment is the criteria we are using to perform unrelated post hoc analysis. In truth there is a degree of circular mathematics in any algorithm used to produce an alignment, but we overcome this by optical verification precisely avoid this.
In this sense the alignment is simply to establish homologous amino acid/nucleotide positions for downstream analysis. Thus optical optimisation, whilst cautioned against in every other area of bioinformatics, is considered useful in alignments because of the risk that a logical error in the alignment is amplified in downstream analysis, potentially leading to an erroneous conclusion.
Algorithms are a problem where indels are concerned and aligning across an indel. It is very easy to make errors here and 'optical analysis' would simple introduce a blockwise deletion to avoid the uncertainty, i.e. if you are not sure just chop it out - better to lose data than risk artefacts on an equivocal judgement call.
I'm not saying don't try, I am saying don't look to overstate the applicability of the metric.


Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ?

The common solution for scRNA-seq is to put cell barcodes and such in read headers and then post-process things with UMItools.

But regarding your actual question, STAR can accept SAM/BAM as input with the --readFilesType SAM PE option (it's SAM for both that format and BAM). You can swap SE for PE if your data is single-end (or is effectively that way, as is the case in many scRNA-seq experiments). Note that in the case of PE experiments that mates must be next to each other. I really don't know what this does with auxiliary tags though, so you would need to see if they're lost before you start processing a lot of data).


Cluster is split in 2-3 locations on tsne plot - Suerat

Your cluster labels come from graph clustering implemented in the FindClusters() function. The resulting clusters are then visualised with a 2D tSNE plot (via RunTSNE() and TSNEPlot()). So, your cluster 13 is not split into three sub-clusters but cells within cluster 13 look somewhat distant from each other on the tSNE plot. Having mentioned distances on a tSNE plot, please see this SE Data Science entry Can closer points be considered more similar in T-SNE visualization? and especially the comment "What I mean is that you can't just use distance in the lower space as a similarity criterion". In your case "the lower space" would be your tSNE plot.

As far as I can see from your code, you are following the widely-adopted "Seurat workflow" and all seem fine including your TSNE plot output. I could not spot the QC step though: apart from the one at object initiation, one usually applies filters for number of genes and % mitochondrial gene contribution (using FilterCells() in Seurat v2 and subset() in Seurat v3).

It would be a good idea to try UMAP instead of tSNE, please see Evaluation of UMAP as an alternative to t-SNE for single-cell data (Becht et al, 2019). Seurat v3 vignettes are now using UMAP and not tSNE any more.

Lastly, please see How to Use t-SNE Effectively to learn about points to take into account when using/interpreting tSNE plots. You will see that you would get surprisingly dissimilar tSNE plots from the same underlying data with different parameters (some are hardcoded in the RunTSNE() function).


Normalizing microarray data for clustering heat map

You see negative values with your function because you're setting the average of each row to 0 and its standard deviation to 1.

In general, I would trust a standard normalization method (rma in this case) more than some random "truncate and then scale the rows" method. Your method isn't even doing any between-array normalization, which is the benefit of rma.


How can i modify FASTA headers in a multi fasta file using BioPython SeqIO

You're very close, but unfortunately for legacy reasons this aspect of the Biopython parser is not very intuitive. The key point is the .id is anything up until the first white space character. One day I hope the .description and .id attributes are replaced by a single attribute like .header. Anyway, here is a solution for your case:

from Bio import SeqIO

to_add = "TEST~~~TEST~~~"

with open("out.fa", "w") as outputs:
    for r in SeqIO.parse("in.fa", "fasta"):
        r.id = (to_add + r.description).replace(" ", "_")
        r.description = r.id
        SeqIO.write(r, outputs, "fasta")

Output on your sample input:

>TEST~~~TEST~~~AAA_a_aa
ATGATG
>TEST~~~TEST~~~BBB_b_bb
GACGAC
>TEST~~~TEST~~~CCC_c_cc
TTGTTG


single stranded cDNA vs double-stranded cDNA

Sequencing also requires single stranded DNA, otherwise the sequencing primers can’t bind. Sure there are steps in the protocol where double stranded material exists, but that’s transient.


I believe you may be able to map your reads, but I don't know how to do that with minimap2.

I recommend running gsnap, which is more SNP tolerant and provides a number of parameters which are likely to help.

For instance, I believe that most aligners will treat 'N' characters as mismatches when you align. GSNAP has parameters to account for this.

--query-unk-mismatch=INT       Whether to count unknown (N) characters in the query as a mismatch
                               (0=no (default), 1=yes)
--genome-unk-mismatch=INT      Whether to count unknown (N) characters in the genome as a mismatch
                               (0=no, 1=yes (default))

It also has a mismatch parameter similar to the one you described for minimap2.

-m, --max-mismatches=FLOAT     Maximum number of mismatches allowed

Try running with the 'genome-unk-mismatch' parameter above (with the masked reference). That may be your best bet. There are also other parameters that may help, but this should be a good start.


What is the R/Bioconductor preferred way of accessing TCGA data?

There isn't a R/Bioconductor person or view, so there isn't a preferred way to access TCGA. Use what suits you.

Bioconductor is build by contributors which are encourage to not overlap in their contributions. Usually the packages have different options/features even if they allow to access to the same data.

In this case, the TCGAbiolinks package is specifically for TCGA dataset while GenomicDataCommons is for "access the NIH / NCI Genomic Data Commons RESTful service". So the functionality they bring to R and Bioconductor is different.

The fact that one project is finishing doesn't mean that the data or the packages build around that project are obsolete. Some functionality different between them is the functions to harmonize and prepare (GDCprepare) data from the TCGAbiolink package, or see a boxplot correlation and AAIC plot to define outliers TCGAanalyze_Preprocessing.


You should never install python packages as root or with sudo.

You're already using conda, so conda install -c bioconda -c conda-forge dendropy.


Parsing MSigDB Supplementary Collections

If you want to separate each set, you can do something like this

dat <- readLines("~/Desktop/scsig.v1.0.metadata.txt")
dat <- gsub("\t"," ",dat)

# get the indeces of each type
standardName <- grep("STANDARD_NAME",dat)
organism <- grep("ORGANISM",dat)
organSystem <- grep("ORGAN_SYSTEM",dat)
pmid <- grep("PMID",dat)

# return the values
standardName <- gsub(pattern = "STANDARD_NAME",replacement = "",dat[standardName])
organism <- gsub(pattern = "ORGANISM",replacement = "",dat[organism])
organSystem <- gsub(pattern = "ORGAN_SYSTEM",replacement = "",dat[organSystem])
pmid <- gsub(pattern = "PMID",replacement = "",dat[pmid])

df <- data.frame(standardName,organism,organSystem,pmid,stringsAsFactors = FALSE)

 head(df)
                                     standardName      organism             organSystem      pmid
1      Fan_Embryonic_CTX_Big_Groups_Cajal_Retzius  Homo sapiens  Central Nervous System  29867213
2  Fan_Embryonic_CTX_Big_Groups_Brain_Endothelial  Homo sapiens  Central Nervous System  29867213
3  Fan_Embryonic_CTX_Big_Groups_Excitatory_Neuron  Homo sapiens  Central Nervous System  29867213
4              Fan_Embryonic_CTX_Big_Groups_Glial  Homo sapiens  Central Nervous System  29867213
5       Fan_Embryonic_CTX_Big_Groups_Brain_Immune  Homo sapiens  Central Nervous System  29867213
6          Fan_Embryonic_CTX_Big_Groups_Microglia  Homo sapiens  Central Nervous System  29867213

This method is very verbose, but it's a start and it can be easily improved with a loop.

 #---- EDIT

Not knowing how the format works, there may be two issues:

1) gene sets with missing values only

STANDARD_NAME Aizarani_Liver_C33_Stellate_cells_2
ORGANISM
ORGAN_SYSTEM
PMID 31292543

2) gene sets with missing value and type

STANDARD_NAME Aizarani_Liver_C33_Stellate_cells_2
PMID 31292543

Case 1 is not an issue, as the code above will simply add an empty cell and go to the next line.
For Case 2 we can use unique to get a list of all define types, create an empty data.frame with those types and have a loop (or sapply) to add the values in the right columns.

We can check the nature of metadata.txt:

 types <- gsub("^(.*)\t.*$","\\1",dat)
 sapply(unique(types), function(t) length(grep(t,types)))

     STANDARD_NAME                ORGANISM            ORGAN_SYSTEM                    PMID       PUBLICATION_TITLE
                257                     257                     257                     257                     257
            AUTHORS                   GEOID            EXACT_SOURCE    EXTERNAL_DETAILS_URL                    CHIP
                257                     257                     257                     257                     257
      CATEGORY_CODE             CONTRIBUTOR         CONTRIBUTOR_ORG       DESCRIPTION_BRIEF RAW_PUBLICATION_MEMBERS
                257                     514                     257                     257                     257

Your are lucky and you don't have to worry about shifting because of missing values.
There is a fixable issue with CONTRIBUTORS, but it is related to the regex, not the file format.


tximport txidcol and geneidcol selecting wrong column

By setting type="salmon" you're overriding your settings for things like txIdCol. The various type options (except none) are presets to set things correctly for the output from various programs like salmon or kallisto. If you need to specify custom column names, set type="none".

As an aside, if I'm interpreting your question correctly you have only two samples in total. You can not perform a meaningful analysis with an unreplicated experiment, so if that is indeed your case then don't waste your time further with this dataset.


Why Do I Get Different Results Each Time From AutoDock Vina?

Vina uses a non-deterministic docking algorithm. See the manual.

The docking algorithm is non-deterministic. Even though with this receptor-ligand pair, the minimum of the scoring function corresponds to the correct conformation, the docking algorithm sometimes fails to find it. Try several times and see for yourself. Note that the probability of failing to find the mininum may be different with a different system.

The manual states that you can set the random seed in order to obtain reproducibility.

I changed something, and now the docking results are different. Why?
Firstly, had you not changed anything, some results could have been different anyway, due to the non-deterministic nature of the search algorithm. Exact reproducibility can be assured by supplying the same random seed to both calculations, but only if all other inputs and parameters are the same as well. Even minor changes to the input can have an effect similar to a new random seed. What does make sense discussing are the statistical properties of the calculations: e.g. "with the new protonation state, Vina is much less likely to find the correct docked conformation".


how to obtain GRCh38 bed and symbol files for HOCOMOCOv11

I guess you're right that you need predictions of transcription factor binding site (TFBS) through genome.
There're two options to get them:

You can generate predictions using any motif scanning tool. E.g. you can use our tool SPRY-SARUS which is able to generate bed-file in necessary format.
In this case, you can choose more or less strict thresholds for each prediction (i.e. more false positives vs more false negatives). Standard thresholds for HOCOMOCO motifs can be found at its downloads page (we, as HOCOMOCO authors, recommend thresholds corresponding to P-value of 0.0005).
I had no experience with diffTF but it seems that authors prepared TFBS predictions and put them on their site. Please, take a look at this page of manual. If you don't need extra-flexibility in the generation of TFBS predictions, it can be the simplest way to proceed.


How can I run a command for multiple files?

Assuming all your vcfs are in the same folder:

for file in *.vcf; do
    grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file}.txt
done


Global alignment between two sequence X and Y with maximum number of identical matches

There are a lot of possible ways to count the number of perfect matches in an alignment, and the details will depend on which alignment tool you are using.
But, if you use the built-in pairwise aligner in BioPython (Bio.Align.PairwiseAligner(), https://biopython.org/docs/1.75/api/Bio.Align.html), then the easiest way is to use the formatted string representation of the alignment, which uses the pipe/vertical-bar character ('|') to indicate perfect matches.
Example code:
from Bio import Align

aligner = Align.PairwiseAligner()
alignments = aligner.align("GAACT", "GAT")
max_matches = 0
for i, a in enumerate(alignments):
    align_str = f'{a}'
    align_str = align_str.split('\n')
    match_str = align_str[1]
    match_count = match_str.count('|')
    mismatch_count = match_str.count('.')
    print(f'Number of matches: {match_count}')
    print(f'Number of mismatches: {mismatch_count}')
    if match_count>max_matches:
        max_matches = match_count
        max_ind = i

best_alignment = alignments[max_ind]

In this example, both alignments have the same number of perfect patches (3), but that may not always be the case, depending on the input sequences and the scoring parameters used. The example just uses the default scoring parameters, but in general, you probably want to choose parameters appropriate for your type of alignment.


The title is misleading as this error doesn't have anything to do with making the font italic. In Seurat v2, FeaturePlot does not return a ggplot2 object by default, so p in your case is NULL. You need to set do.return to TRUE in the FeaturePlot call.

You should instead do:

library(Seurat)

p <- FeaturePlot(pbmc_small, head(pbmc_small@var.genes), do.return = TRUE)

for(i in 1:length(p)) {
  p[[i]] <- p[[i]] + theme(plot.title = element_text(face = 'italic'))
}

cowplot::plot_grid(plotlist = p)


Parallel in PLINK for linear association for SNP effects

Try the latest plink 2.0 alpha build for this.  It’s an incomplete program, but its implementation of —linear is far better than v1.9’s.


It might be used differently in different contexts, but generally speaking, in my world - allelic imbalance is when there's a difference in the level of gene expression from different alleles, usually through genetic (e.g. a variant in a promoter) or epigenetic mechanisms (e.g. one copy silenced, as in imprinted regions).

There are a few blog posts here with some more information, and plenty of papers (google "Allelic imbalance") if you want to dig deeper.


How does MinKNOW classify 1D reads as "pass" or "fail"?

The classification is indeed based purely on comparing a mean to a threshold... but it's not the mean of the Phred scores that represent base quality on a logarithmic scale, but rather the mean of the underlying quality values prior to being transformed into Phreds.

If I compute the "mean" quality for a read in this way - by converting each Phred in the read back to a probability, taking the mean of those probabilities, and then converting that mean probability to a Phred - I find that whether a read is a pass or fail is determined by comparing that score to a threshold of 7:

A mean quality score computed in this way is also included in the mean_qscore_template column of the sequencing summary file. If you open that file in a spreadsheet program (it's just a TSV) and sort by mean_qscore_template, you will find that the passes_filter column is always FALSE for mean_qscore_template values below 7 and always TRUE for values above 7.

(Hat tip to Wouter De Coster for successfully guessing that this was the answer in the comments and pointing me on the right path.)


How to query the Human Microbiome Project (HMP) to find all subjects with both 16s and WGS workups?

After trying a few different queries, I don't believe there is a more specific one than this and I found that by searching and clicking through all of the results given by this query that eventually I found patients that fit the criteria of having both WGS and 16s.


Transcript Coordinate Ranges to Genomic Coordinates

This task has the same "flavor" as many I've done before, but each case is so subtly different that it's impossible to write a generalized tool that will work correctly for all circumstances.

R isn't my wheelhouse, but I was able to throw something together pretty quickly with standard Python. This should get you at least 95% of the way to where you're going. :-)

#!/usr/bin/env python3

from argparse import ArgumentParser, FileType
from collections import defaultdict
from re import search, sub

cli = ArgumentParser()
cli.add_argument('txspace', type=FileType('r'))
cli.add_argument('chrspace', type=FileType('r'))
args = cli.parse_args()

exons = defaultdict(list)
for line in args.txspace:
    if line.startswith('#') or line.strip() == '':
        continue
    values = line.split('\t')
    featuretype = values[2]
    if featuretype not in ('CDS', 'exon'):
        continue
    txid = values[0]
    exons[txid].append(values)

for line in args.chrspace:
    values = line.split('\t')
    if len(values) != 9:
        print(line, end='')
        continue
    featuretype = values[2]
    if featuretype not in ('gene', 'mRNA'):
        continue
    print(line, end='')
    if featuretype == 'mRNA':
        txid = search('ID=([^;\n]+)', values[8]).group(1)
        start = int(values[3])
        txexons = exons[txid]
        for exon in txexons:
            exon[0] = values[0]
            exon[3] = str(int(exon[3]) + start)
            exon[4] = str(int(exon[4]) + start)
            exon[8] = sub(r'Parent=[^;\n]+', 'Parent=' + txid, exon[8])
            print('\t'.join(exon), end='')


I'd recommend using R and either merge or one of the joins. Assuming your data is not too large, data.frame + read.delim should suffice to hold the two datasets. If not, try data.table + fread..


Visualize network (graph) using "line" layout

That plot is called an Arc diagram
you can google that and use the package or program that best suits you.

Here are 2 examples of the use of Arc diagrams in R:
Example 1:

Example 2:


Bedtools wrongly indicates zero coverage

If this is your output :

chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000

Then I guess your bedfile looks like this (the first 6 fields) :

chrom : chr1
FeatureStart : 12595
FeatureEnd : 12802
name : ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242
score : .
strand : ref|DDX11L1

The issue I see here is that you are giving the wrong information to the strand field of your bedfile. It should be either +, either - or . if you don't want to consider the strand information.

Hope it helps !


What do the symbols mean in minimap2's gap cost equation?

Here is the translation:

$l$ = gap length
$q$ = initial gap penalty for shorter gaps
$e$ = short gap extension penalty per every additional gap
$\tilde{q}$ = initial gap penalty for longer gaps
$\tilde{e}$ = long gap extension penalty per every additional gap

To understand the first of the two expressions being minimized, see the Wikipedia section on affine gap penalties.

Notice the similar notation structure

$$A + B \cdot L$$

where

$A$ is the known gap penalty,
$B$ is the gap extension penalty, and
$L$ is the length of the gap

To understand the other variables with the tildes on top, see the text below Equation 2 (emphasis added to key explanation):

  On the condition that $e > \tilde{e}$,
  it applies cost
  $q + |l| \cdot e$
  to gaps shorter than
  $[(\tilde{q} - q)/(e - \tilde{e})]$
  and applies
  $\tilde{q} + |l| \cdot \tilde{e}$
  to longer gaps. This scheme helps to recover longer insertions and deletions (INDELs).

In other words, we make the assertion that we should penalize longer gaps differently than shorter gaps.

As a bonus, here's a visual understanding of why we use the minimization function on those two expressions using Python and some made-up (and probably not biologically useful) penalty terms. Note the requirement that $e > \tilde{e}$.

# Import libraries
import numpy as np
import matplotlib.pyplot as plt

# Setup
l = np.array(range(0, 15, 1))
y_short = eval('10 + 5*l') # Q=10, E=5
y_long = eval('20 + 3*l')  # Q=20, E=3

# Plot
plt.plot(l, y_short)
plt.plot(l, y_long)
plt.legend(["Short Gap Penalty", "Long Gap Penalty"])

Using our example and the expression $[(\tilde{q} - q)/(e - \tilde{e})]$, it will apply the first expression in the minimization function (the one without the tildes) when the gap length is shorter than 5.

$$
[(\tilde{q} - q)/(e - \tilde{e})]
=
[20 - 10)/(5 - 3)]
=
10 / 2
= 5
$$

We can see from our plot that the two lines intersect at 5, which we've verified through some arithmetic.


Whenever you have batch or other nuisance variables driving variation you're going to get odd networks. You tried using SVA, which is a good method to use, but it seems in this case that it's able to optimally handle this dataset. I also strongly suggest that you use pickSoftThreshold() to choose a softpower, since that may produce better results (though since you already have negative values I'm not holding by breath). In general, you're going to need to play with either SVA or filtering to remove the batch-related genes that are causing problems.


Is it important to filter out poor quality cells before performing an integration analysis on single cell RNA sequencing data?

It is absolutely necessary to remove low quality cells: In the case of CCA (and this applies to other "integration" or "data alignment" methods as well), one would need to use "anchors", basically same type/state/kind of cells from the samples to be "integrated" and are used to "align" the different samples. In a scenario where you have a set of cells whose distinguished feature is high mitochondrial content due to apoptosis irrespective of their origin (epithelial, T cell, ...) and such cells are picked up as so called anchors, the whole integration effort will be off. Even when not using integration methods, bad quality cells quite often form their own clusters (highly skewed biology) not contributing anything to analyses and are often removed at the downstream steps anyway.

And when it comes to published work, unfortunately it comes down to the backgrounds of the reviewers. I tend to think that there are way too many single cell studies than experts with required bioinformatics and statistics so that a thorough review can be conducted. And lack of a paragraph on QC does not mean no QC has been performed anyway.


How to filter Ensembl cDNA and ncRNA FASTA files by primary assembly?

The simplest method is to download the GTF file for GRCm38 and filter that. You can then use one of the many tools out there (bedtools getfasta (make sure to enforce strandedness and possibly use the --split option, if you happen to use a BED12 file), gtf_to_fasta from tophat2 (it's the only reason to keep tophat2 installed), etc.) to make a transcriptomic fasta file.

Alternatively, Ensembl BioMart can be used to directly download a transcriptomic fasta file using a filter to select only those transcripts on pre-defined chromosomes. Other filters can be applied as well. If you don't already have a GTF file downloaded and one of the previously mentioned packages installed, then this may be the simplest way to get what you want.


Searching for start and stop codons for protein sequencing of contigs

If your string is 'GGTAGCCAATTATGAACGAATAGGAC' then I would assume that you could separate this into coding and non-coding regions such as:
GG*TAG*CCAATT*    (non-coding)
ATG*AACGAA*TAG    (coding region)
*GAC              (non-coding)

It starts with a start codon ATG and ends with an end codon TAG.
Your first event position 12 does not make sense to me. What does * mean? It seems of that your sequence starts with GG which is not a triplet? Perhaps did you mean GGG?
I'm the author of node based DNA visualisation AminoSee. It has a bug in that it only renders half right: because DNA can be read backwards, and my program only does it forwards. I actually came to BioInformatics today to try and get my head around reading backwards and how that might work.
AminoSee only transcodes 3 base pairs at a time, and discards triplets that aren't, well, 3 base pairs longs. I'm self taught and so I guess unaware if it's possible to have any sequence that is not a multiple of 3. Perhaps it's possible. I assumed not. My program is more of an art piece that a science tool, but I plan to use it one day to perform a "visual diff" between two species (Humans and Chimp). See AminoSee  DNA Visualisation over at aminosee.funk.nz


TRAL 0.3.5 has bugs running and parsing Phobos. You can follow progress towards a fix on this issue: https://github.com/acg-team/tral/issues/16

Until the issue is fixed, you should remove Phobos from the repeat detector list in your tral config file (~/.tral/config.ini) and use a different repeat detector.

[sequence]
    [[repeat_detection]]
        DNA = TRED, T-REKS, TRF, XSTREAM


How to retrieve the subcellular location info from uniprot?

Since you are already on UniProtKB, the simplest way will be if you add this information to the output table and download it as a that you can import and filter in R.

To get additional annotations, click the button labelled "Columns" on top of the UniPortKB website.

This allows you to add a wide variety of annotations such as subcellular location to the output. Tick the checkmark for subcellular location and save on the top right.

Then you can download the proteome-wide annotation for example as a tab-separated list using the "Download" dialog next to the "Columns" selection. This list can be imported into R.


Dataset: Locations of regulatory sequences in the human genome?

This can be done quite easily using Ensebl's BioMart.

Choose the Ensembl Regulation database:

Select the "Human Regulatory Features" dataset:

That's basically it right there, just click on "Results":

Export to file and click "Go":

This will download a file called mart_export.txt which looks like this (I chose TSV for tab separated values):

$ head mart_export.txt
Chromosome/scaffold name    Start (bp)  End (bp)    Feature type
18  76429380    76430144    Open chromatin
8   66405962    66406502    Open chromatin
4   61184401    61184600    CTCF Binding Site
X   40733600    40737000    Promoter
5   97407001    97407200    CTCF Binding Site
X   73946201    73946600    Promoter Flanking Region
15  19948201    19949200    CTCF Binding Site
5   11302601    11303143    Open chromatin
2   208407801   208408000   CTCF Binding Site

As you can see the format is pretty straightforward. The first column has the chromosome or scaffold name (it should always be a chromosome for human), the second and third are the star and end positions and the last field is the type of region.

You can also choose to limit by region type by choosing a filter (click the "Filters" link) in BioMart before downloading or by simply parsing the file once you've downloaded it.


How to assess if this strategy is good to have some phylogenetic information?

The method you are proposing is called a distance method, because you are constructing a distance matrix and then using a clustering algorithm to make a tree. The things to watch out for are:

Multiple substitutions, traditionally this is take care of using a Jukes-Cantor correction. You could argue this is less of a concern because this is a 64 substition matrix model rather than just 4 or 20 substition matrix for nucleotides or amino acids.
Use the neighbour-joining clustering algorithm. Do not use a standard hierarchical clustering algorithm such as UPGMA. I'm sure R must have the nj in something like R's ape, pretty sure scipy will have something like that
Code for bootstrap resampling. That is a method for measuring robustness

The codon model for phylogenetic trees has been implemented in MrBayes, I suspect it will be in phyloBayes. This method will be Bayesian and you should run this is parallel with your own solution.
Concerns
Codon models are very prone to mutation saturation (which is complicated) so you must only use it for closely related sequences. If e.g. 3rd codon distances start exceeding 70% codon models go kaput, because it will select a different codon, which is a completely separate part of the matrix - completely by chance. So whilst mutation saturation in a nucleotide model (4 substitution matrix) would invalidate the tree in theory, they would be sufficient 'good' signal for key phylogenetic signals to be identified, this wouldn't happen for a codon model.
I recognise the 'concerns' might sound a bit confusing at present, but just to reiterate, codon models need sequence data that is really closely related. For this data they excel and give excellent bootstrap support.


Depends on the software to generate the segmentation file;

For example I have used scatngs for segmentation that for sample returns

Segment number
Chromosome
Start position (origin-1)
End position (origin-1)
Major copy number—normal
Minor copy number—normal
Major copy number—tumor
Minor copy number—tumor

By  log2(Major copy number—tumor +  Minor copy number—tumor /ploidy of each sample) you would get segment mean

For number of probes

The number of probes depends on where the data came from.
1) From SNP arrays:  Number of array probes within each segment.
2) From exon capture data:   Probes are the ends of exons and so the number of probes are equal to the number of exons  X 2.
3) From whole genome sequencing:  Probes are the boundaries of bins used to generate the copy number data and so the number of probes are equal to the number of bins +

For example in my case for scatngs

Each SNP is a marker, and so you have to count how many of these SNPs are on each segment

I have asked this question in

Extracting the number SNP in each range


Read length distribution from FASTA file

If you want something quick and dirty you could rapidly index the FASTA with samtools faidx and then put the lengths column through R (other languages are available) on the command line.

samtools faidx $fasta
cut -f2 $fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)'

This outputs a statistical summary, and creates a PDF in the current directory called Rplots.pdf, containing a histogram.


Intersecting a data frame based on the multiple columns

reformat df2 so that you have chr1 is instead of 1, then you basically have two bed files.  Then use bedtools intersect.  https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html


How do I convert an Excel file (.xlsx) to an .smi file format?

So a smi file is simply a text file formatted SMILES tab name newline.
Therefore what you need to do is

save for safety
make a new sheet
copy the SMILES column and the name column
save as tab separated file
in Finder/Nautilus/File Explorer rename the file extension to .smi


How to deal with heterozygosity during polishing of genome assembly based on long reads?

A few possibilities:

Falcon

Try falcon and falcon-unzip. These are designed exactly for your problem and your data: https://github.com/PacificBiosciences/FALCON

Not Falcon

If you think you have assembled haplotypes (which seems reasonable to expect given enough coverage), you should be able to see the two haplotypes by just doing all pairwise alignments of your contigs. Haplotypes should show up as pairs of contigs that are MUCH more similar (even with a lot of between-haplotype divergence) than other pairs. Once you have all such pairs, you can simply select one of each pair to polish.


Mathematical notation for formulating a rank score

Here is my solution to the problem. I am posting it here in case anyone else come up with the same idea but did not know how to formulate it in mathematical notation!

  $$ \psi_i = \sum_{i=1}^9 S_i $$

  where $S_i$ is the score function and defined as:

  $$
S_i = \cases{
  0 & \text{if $\theta(x_i) < Q_1$ or $\theta(x_i) = $ benign/neutral} \\
  0.5 & \text{if $Q_1 \leq \theta(x_i) < Q_3$ or $\theta(x_i) = $ possibly damaging/uncertain} \\
  1 & \text{if $\theta(x_i) \geq Q_3$ or $\theta(x_i) = $ damaging}
}
$$

  The $\theta(x_i)$ is the pathogenicity or conservation score for variant $x$ as defined by model $i$ and $Q$ denotes the quartile range for scores from M-CAP, CADD, GERP and Phylop models.


Question of Padj value (very high, infinite) got from dds

The fold changes for those are huge (on the order of 32x for the largest), so it's unsurprising for the adjusted p-values to be absurdly small.
There's no reason to treat them differently.


How to check if a given gene is expressed in a group of microarray samples if I do not have control group to compare with?

You ask about which "genes are expressed" and then you mention "if a gene is up or down regulated". These are different, and given your application I think what you actually want to know if these marker genes are expressed. This is not a question of up/down-regulation relative to a control, so you do not need control samples for this.
It seems most of your datasets are Affymetrix, and there are R/Bioconductor packages that will give you the probability a gene is expressed. One such I found is frma (https://bioconductor.org/packages/release/bioc/html/frma.html), which returns an object "weight" that is "a vector of weights which roughly correspond to the probability of expression for each gene." I haven't worked with microarrays in a while, but I recall many other R/Bioconductor packages offering this functionality.
A detection p-value is calculated for Illumina array probes in Illumina's native BeadStudio software.


Split FASTQ and matching BAM into matching chunks

Hmm, it's hard to think of a super efficient way of doing this (assuming the files aren't ordered the same - if they are then this whole answer is basically redundant). And also assuming the read ids for both files aren't a perfect intersection.

Off the top of my head you probably want to build a set of read ids for the fastq file and another for the bam. Some python code to get started on that:

import pysam
import itertools

def get_read_id_fastq(ref_path):
    """Extracts the read ids from a fastq file."""
    read_ids = set()
    with open(ref_path, 'r') as ref:
        for line in ref:
            if line.startswith('@'):  # i.e if line is header
                # split the line on spaces, take the first element, remove @
                read_id = line.split(' ')[0].replace('@', '')
                read_ids.add(read_id)

    return read_ids

def get_read_id_bam(ref_path):
    """Extract the read ids from a BAM file."""

    read_ids = set()
    with pysam.AlignmentFile(ref_path, 'r', ignore_truncation=True) as ref:
        for read in ref:
            # query_name is the query template name
            read_ids.add(read.query_name)

    return read_ids

fastq_ids = get_read_id_fastq(fastq_path)
bam_ids = get_read_id_bam(bam_path)

Then take the intersection of these two sets and that's your common read ids.

common_ids = fastq_ids & bam_ids

The next part will be a bit more involved. You will have to iterate through each file one at a time. I would suggest that for the first one you iterate through, create a running dictionary with the read id that was written as the key and the 'chuck number' it was written to as the value. You could create a cycle for your chunk size to manage this effectively

chunk_cycle = itertools.chain(*zip(range(chunk_size)*len(common_ids)))
write_idx = {}

The next part will probably require you to put in some hard-fought times getting it to work. I'll put in some rough pseudo-code to give you an idea.

for line in file:
    read_id = # get line read_id
    if read_id in common_ids:
        chunk_num = chunk_cycle.next()
        write_idx[read_id] = chunk_num
        file_to_write_to = 'out.{}.bam'.format(chunk_num)
        # open this file or write to it if already open

When you go to do the next file when you find a read id is in the common set, you look up it's value in write_idx and this will give you the chunk number to write to.

The reads wont be in the exact same order between the two files, but you could sort later if you needed it (not sure it would matter?).

Hope this helps. Sorry I couldn't give you more, but this should give you a head start hopefully.


How to quickly and robustly convert between mmCIF and PDB?

Converting a PDB file to some version of an mmCIF file is always possible provided the PDB file doesn't have any "problems". However grouping identical molecules into the same entity is non-trivial so the Biopython and BioJulia mmCIF writers (which I wrote) do not try and do this, simply assigning new entites for each molecule. These files will also be missing other things like the metadata parts of the mmCIF dictionary not involving the atoms.

In Biopython:

from Bio.PDB import PDBParser, MMCIFIO
p = PDBParser()
struc = p.get_structure("", "file.pdb")
io = MMCIFIO()
io.set_structure(struc)
io.save("file.cif")

In BioJulia:

using BioStructures
struc = read("file.pdb", PDB)
writemmcif("file.cif", struc)

Going the other way is, however, not always possible. You can run out of atom serials or chain IDs to use in the output PDB file (this is one reason for the switch to mmCIF files in general). The gist you have linked is a way to get round this.

If this is not a problem for your structure you can invert the above snippets to get a decent PDB file - see the relevant docs. Note that in both cases it will use the auth_seq_id and auth_asym_id fields to get the residue number and chain ID rather than the relabelled label versions. The test suites of Biopython and BioJulia do test these conversions pretty thoroughly.

There is also pdb-tools if you do not want to install Biopython, though both involve a simple pip install so it's not much difference.

For some speed benchmarks of various software see here and here.


Can I model technical replicates in DESeq2?

Practically speaking, there's no way to include the technical replicates in that design (in DESeq2 at least). Your concern regarding inflating the power is exactly correct and the only way to combat that would be to add a pairing factor like one might do with case-control or tumor-normal studies. That is, something like:

  group libraryPrep sample
1    WT           A      1
2    WT           B      1
3   MUT           A      2
4   MUT           B      2
5    WT           A      3
6    WT           B      3
7   MUT           A      4
8   MUT           B      4

Here sample pairs the technical replicates together. However this ends up either being rank deficient or in the end not more informative.


The state, limitations and comparisons of large variant stores

An epic question. Unfortunately, the short answer is: no, there are no widely used solutions.

For several thousand samples, BCF2, the binary representation of VCF, should work well. I don't see the need of new tools at this scale. For a larger sample size, ExAC people are using spark-based hail. It keeps all per-sample annotations (like GL, GQ and DP) in addition to genotypes. Hail is at least something heavily used in practice, although mostly by a few groups so far.

A simpler problem is to store genotypes only. This is sufficient to the majority of end users. There are better approaches to store and query genotypes. GQT, developed by the Gemini team, enables fast query of samples. It allows you to quickly pull samples under certain genotype configurations. As I remember, GQT is orders of magnitude faster than google genomics API to do PCA. Another tool is BGT. It produces a much smaller file and provides fast and convenient queries over sites. Its paper talks about ~32k whole-genome samples. I am in the camp who believe specialized binary formats like GQT and BGT are faster than solutions built on top of generic databases. I would encourage you to have a look if you only want to query genotypes.

Intel's GenomicDB approaches the problem in a different angle. It does not actually keep a "squared" multi-sample VCF internally. It instead keeps per-sample genotypes/annotations and generates merged VCF on the fly (this is my understanding, which could be wrong). I don't have first-hand experience with GenomicDB, but I think something in this line should be the ultimate solution in the era of 1M samples. I know GATK4 is using it at some step.

As to others in your list, Gemini might not scale that well, I guess. It is partly the reason why they work on GQT. Last time I checked, BigQuery did not query individual genotypes. It only queries over site statistics. Google genomics APIs access individual genotypes, but I doubt it can be performant. Adam is worth trying. I have not tried, though.


How to create a list of differentially expressed (DE) genes after normalization with RUVSeq?

You do the normalization before running your edgeR. The purpose of RUVg is to remove "Remove Unwanted Variation Using Control Genes". In your code, you ran edgeR and then normalize the data using RUVg, which is only going to return you the normalized counts.
Using the example dataset in vignette:
library(RUVSeq)
library(zebrafishRNASeq)
data(zfGenes)
filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]
genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]
spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]

x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered),
                           phenoData = data.frame(x, row.names=colnames(filtered)))
set <- betweenLaneNormalization(set, which="upper")

set1 <- RUVg(set, spikes, k=1)

You can look at it, it's an expression set with counts etc, not results:
set1
SeqExpressionSet (storageMode: lockedEnvironment)
assayData: 20865 features, 6 samples
  element names: counts, normalizedCounts, offset
protocolData: none
phenoData
  sampleNames: Ctl1 Ctl3 ... Trt13 (6 total)
  varLabels: x W_1
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:

You run edgeR now on the results of RUVg:
design <- model.matrix(~x + W_1, data=pData(set1))
y <- DGEList(counts=counts(set1), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
topTags(lrt)


Yes, Prokka will consider all six-frames, the gene prediction is actually done with Prodigal, paper:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2848648/

github:
https://github.com/hyattpd/Prodigal


Can I create a CRAM file with a relative reference path?

In your CRAM file, make all of the reference sequences present in the ./ directory. Then make a bash script to make relative links to your reference for whatever mount you are using at the time - ie ln -s TARGET LINK_NAME.

Now your CRAM is portable and you just have to tailor the local working dir to it.


replacing SNPs with missing calls with a specific string

A possible approach in plain python:

#!/usr/bin/env python3

import sys

def both_zero(freq1, freq2):
    return freq1 == "0" and freq2 == "0"

with open(sys.argv[1], "r") as snp_file:
    for line1 in snp_file:
        # Consume one extra line to get the frequencies
        # for the alternative allele
        line2 = snp_file.readline()
        # "zip" the frequencies in pairs
        freq_pairs = zip(line1.strip().split(), line2.strip().split())
        recoded_pairs = [
            ("9", "9") if both_zero(freq1, freq2) else (freq1, freq2)
            for (freq1, freq2) in freq_pairs]
        # "unzip" the recoded frequencies into separate objects
        # (see https://stackoverflow.com/a/19343/1878788)
        (recoded1, recoded2) = zip(*recoded_pairs)
        print(*recoded1, sep="\t")
        print(*recoded2, sep="\t")

Calling it with your example snp frequencies file as command-line argument:

$ ./recode_absent.py snps.txt
0.65    0.722222        1.0     0.75    9
0.35    0.277778        0.0     0.25    9
9       0.666667        0.75    0.5     0.5625
9       0.333333        0.25    0.5     0.4375


What's your reference for the definition of STRs? I think it is still ambiguous among the community. Wikipedia states motifs of >=2 bp. However other references include homopolymers as well:
Source 1, source 2 and source 3.


How to assess the quality of assembled .fasta genome files?

To expand on the comment by Devon Ryan, what you really need are (A.) bam alignment files for each individual and (B.) VCF files for each individual showing all of their genomic variants (with respect to the reference human genome).  You mentioned that you had fastq files and that you converted them into fasta files.  Do you know if the bam files and VCF files are also available for download anywhere?

If the bam files are not available you will have generate them using sequence aligner software (for example: bwa-mem).  With the bam files you can calculate the average depth of coverage (and many other things).

If the VCF files are not available then you will have to generate them from the bam files using variant calling software (for example, the GATK package from the Broad Institute).

In terms of assessing the quality of the sequence runs, a commonly applied tool is name FastQC, but once again, it uses fastq formatted files as input, not fasta files.

The main use for fasta formatted files these days seems to be for compact storage of various reference genomes from different genome sequencing projects (using gzip or similar type of compression).

With the VCF files for each individual in hand you should be able to begin addressing your other questions regarding the frequency of shared variants, etc.


What are all the reference files produced by bwa index, and are these dependent upon whether the reference is zipped?

You'll get the exact same index (the amb, ann, bwt, pac and sa files) whether the reference is gzipped or not. BWA also makes its own packed reference sequence (the .pac file) so you don't even need the genome around after you index.


Best way to "fish" long reads that match a query sequence

I have seen the blasr aligner typically used for aligning PacBio reads to a reference.  Is there a reason why blasr would not be a good choice? I guess the very short sequence you are trying to align to might be an issue.

Minimap2 is capable of mapping pacbio reads, and can also do pairwise mapping to see if reads overlap with each other, so this function might help you:
From https://github.com/lh3/minimap2

minimap2 -x ava-pb  reads.fq query_seq.fq > ovlp.paf