The * character has a special meaning in regular expressions, it means "0 or more". So when you write {6,}*, the regex engine tries to apply the "0 or more" to the {6,*} and gets confused:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.6/re.py", line 222, in findall
    return _compile(pattern, flags).findall(string)
  File "/usr/lib/python3.6/re.py", line 301, in _compile
    p = sre_compile.compile(pattern, flags)
  File "/usr/lib/python3.6/sre_compile.py", line 562, in compile
    p = sre_parse.parse(p, flags)
  File "/usr/lib/python3.6/sre_parse.py", line 855, in parse
    p = _parse_sub(source, pattern, flags & SRE_FLAG_VERBOSE, 0)
  File "/usr/lib/python3.6/sre_parse.py", line 416, in _parse_sub
    not nested and not items))
  File "/usr/lib/python3.6/sre_parse.py", line 619, in _parse
    source.tell() - here + len(this))
sre_constants.error: multiple repeat at position 10

In order to use special characters literally, as simple strings, you need to escape them by adding a \ before them. So try:

import re
myprotein = ">Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD"
protein = re.findall('M[A-Z]{6,}\*',myprotein)
print(protein)
['MVTGSCQTLVQDKV*', 'MIPPCVFHPIYPLFTWSPLHPIISGCF*']

To get the output you asked for while reading from a file, you can do:

import sys
import re

proteinFile = sys.argv[1]

with open(proteinFile, 'r') as f:
    for line in f:
        ## Split the long input line into the list 'sequences'
        sequences = line.rstrip('\n').split('>')

        for i, k in list(enumerate(sequences)):
            ## Skip the first element which will be empty.
            ## there must be a pythonic way of doing this, but I'm
            ## a Perl hacker, so I don't know it.
            if len(k) == 0:
                continue

            seqName, seq = k.split(' ')
            ## Get the frame
            frame = seq[0]
            ## Remove the frame from the sequence
            seq = seq[1:]
            # cleave sequence
            peptides = re.findall(r'(M[A-Z]{6,}\*)', seq)

            for idx, peptide in enumerate(peptides):
                print(">%s.%s.%d\n%s" % (seqName, frame, idx+1, peptide) )


Regarding canonical analysis of principal coordinates (CAP) for linking microbial data with elements

I don't think that you can say that the 25% can be said to relate to the chemical variable specifically; rather, that 25% is the variance that is explained by all constrained axes. While these axes are drawn in such a way that they maximize variance associated with your independent variables, I do not think that they are directly interpretable as being 100% due to the elements.
The number of constrained axes is a parameter chosen as an input to the CAPscale method:

An extremely important point is to determine how many PCO axes should be retained; i.e., what should be the choice for $m$? One might be concerned that the method could suffer from a poor choice for $m$. If $m$ is too small, then there may be some ecologically important information in the data that will not be included in the canonical analysis. On the other hand, if $m$ is too large relative to $N$, then a misleading canonical plot could result. For example, we cannot set $m = N - 1$, for then the squared canonical correlations would all be equal to 1.0, because we have $N$ points!

So the value of the variance there will vary depending on how many constrained axes you use. If I understand right, you could obtain 100% distance variance explained by choosing as many constrained axes as you have observations.
I think that it is worth looking into how others have solved this problem. For example, one has to first decide whether or not this 25% is an unexpectedly large number, i.e. compute  p-values (I am not crazy about p-values in most contexts but here they could serve as a helpful heuristic), as suggested by the CAPscale paper:

As a further option, one may test the hypothesis of either (1) no significant differences in multivariate location among groups (CDA), or (2) no significant relationship with quantitative environmental or other variables (CCorA). This is done by using the trace statistic (sum of canonical eigenvalues (sum of squared canonical correlations; see the Appendix) and obtaining a P value by permutation (e.g., Anderson 2001b).

This seems to be directly implemented as the anova.cca() function in the vegan package (I believe that they have decided that CAP is a flavor of CCA for that package). I understand the motivation to get % variance as it's intuitive, but in complex statistical analysis of distances it is a little bit less clear what it really means. Some kind of comparison to a null hypothesis is probably necessary to say whether it is meaningful.
For how others have approached similar questions you can see this discussion. For other implementations that may be more user-friendly, you might consider the vegan package, which is a richer toolbox for ordination in ecology (even if it is less focused on metagenomics). Here is the CAPscale function in vegan. The phyloseq authors have correctly decided to build large pieces of phyloseq out of vegan, so I think that you can pass phyloseq back into vegan functions and e.g. use those analyses directly on your cap_ord object. But see that link for more info.
For a little more easy-to-understand background on statistics of distance methods (particularly in microbiomes), you might see here.


how to generate structure file from blast results of an input sequence

If you the DNA/RNA crystal structure:

Simply change the nt in NCBIWWW.qblast("blastn", "nt", sequence_data) to pdbnt to get the nucleotide sequence in a chain.
The PDB code will be the accession ID with an underscore plus chain letter.
Do note that the PDB structure may be offset.

If you want the crystal structures of the encoded protein:

blastx (translation) on the pdbaa database. As above the PDB code will be as above.

If you want a model of the encoded protein:

Check Swissmodel for a pre-existing close threaded model, else check out I-Tasser (best), Phyre (fast).

If you want a model of the DNA/RNA:

You are running blast without special settings —its not a short motif—, so these things may be long, so do be weary of any result. There are a few RNA folding prediction tools online (e.g. RNAFold).


Make ipyrad use cuda-enabled NVidia card on Ubuntu

You would have to get CUDA Python set up on your system and likely modify ipyrad to actually use it (for example, by incorporating numba). So while this is possible, it will require significant effort on your part. In general, there are very few bioinformatics tools that support GPUs.


Interpreting 0x200 flag in bwa-mem alignments

bwa mem is using exactly the same meaning of the 0x200 flag as every other program, including picard. Don't blindly assume that that entry in the header file relates to the output flags in the SAM output. If you search the code, you'll see that bwa is using this internally as a flag for "this alignment is soft-clipped", which it handles differently for primary and secondary alignments (hard-clipping the latter).

For what it's worth, bwa mem (at least the version I just tried) doesn't set the QC failed flag. This behavour is actually pretty common among aligners, since Illumina has changed how they represent QC failure at least once.


Non-exome variants called from whole-exome sequencing data

These are called off-targets due to unspecific binding of DNA to the capturing beads during library preparation. This is normal and expected since capturing is not a perfect process. It is correct to align against the whole genome but you should simply limit variant calling to exons (or basically the on-target regions) and ignore everything else.


You're right about trying to get the R environments to be as similar as possible, as easily as possible. It's not always necessary to have the exact environment to get the same results, but it definitely helps if you find yourself struggling to replicate results.

General tips for ensuring reproducibility, not specific to R

See Sandve et al., PLoS Comp Bio, 2013. They mention a number of things that you do, like recording seed values and session info. I'll highlight a few points:

Environment: Include a small script that lists all the necessary dependencies

This should include all version numbers of languages and packages used
Better implementation: make it a script that automatically installs dependencies

Data: Share raw and/or processed data where possible

At the very least, details on how to obtain the data you used
This can be an accession number in a database, or a link to the source material
Better implementation: again, make it a script that automatically downloads the raw data you used

Code: Comment and document your code

People (yourself included) need to understand what you've done and why
Additionally, your work may be reproducible, but it may not be correct. If you (or anyone else) attempts to fix a problem in your code, they will have a difficult time doing so if they can't understand what your code is doing

Code: Use version control

Services like GitHub and BitBucket make it simple to share your exact work with collaborators
Version control also ensures that you have a history of all the changes done to your work

No manual manipulation of data

If something is done to your data, write it down

R-specific tips

The packrat package, as you mentioned in your answer
Using a number of custom functions for analyzing your data? Compile those functions into an R package

Allows others to easily load/unload the exact functions you used through tools like CRAN, BioConductor, or devtools::install_github()
Makes it easy to see whether these functions are working as intended, since they are separated from other parts of the analysis

Make unit tests for your work with the testthat package

Helps ensure that your code is working as intended, and gives others confidence in your work as well
Write the tests using testthat::test_that(), and check that they pass using R CMD check or devtools::check()
This works especially well with R packages
This also works especially well with TravisCI, since you can automatically have your code get tested after you change it and make a commit

Consider Docker images for your environment

This bundles everything about the R environment together into one easily executable image (and it's not just limited to R)
These images can be easily shared through Docker Hub
Note this might not always be the best solution, since few people know Docker and it's not always simple to use. Docker also requires admin privileges, which may not be possible for work on computer clusters
There is lots of current discussion about Docker in science. See Boettiger, 2014, for an example, using R


What is and how to detect a dephased read

Phasing in Illumina data means that the read you get is shifted one base before or after what the sequence really should be.  They must be expecting some constant 6 bases at the end of each read, and they are throwing out all reads that don't have those 6 bases.  This would include N's, but a small deletion or insertion would also cause the last 6 bases to not match.  This is in theory detectable and fixable; the dropSeq protocol for instance catches instances where improperly short cell barcodes shift the UMIs, but these authors look like they decided not to try and bother, and just throw out everything that didn't match their exact expectations.


shapeit4 is currently (probably) the most accurate and quickest phasing algorithm which can account for sequencing reads in the model. They suggest you extract the phase information from the reads using WhatsHap first however. The guide for WhatsHap can be found here and has some info on using data processed using GATK.


gene-level versus transcript-level analysis

There are a variety of reasons people use gene-level quantitations.

Transcript-level differences are difficult to biologically interpret. Let's be honest, few groups are likely to put in the work required to determine what these might mean. Most genes have at least some level of characterization, so it's much more biologically tractable to think in terms of gene-level changes. Granted, "genes" are largely artificial constructs and in an ideal world one would always think in terms of transcripts, since they're what's really changing. But we don't live in an ideal world.
Statistical power increases with counts, so if you're splitting your counts among a handful of isoforms you could be missing something. In practice this isn't a big problem, since there are typically obvious major isoforms, but that won't always be the case.
"Mo transcripts mo problems", to channel the bioinformatics notorious B.I.G. It's common practice for people to use a subset of transcripts when using tools like Kallisto and Salmon, since the results usually end up being better. But which subset of transcripts should you use? Those in refseq? Gencode basic? A subset of one of those? Do we really understand all the consequences of choosing one over the others? (No, we don't) We all have our own preferences, but this isn't a solved problem.


Disclaimer, I am not an expert on admixture, answer is provided without any warranty.

It depends what you want to figure out. Why do you want to categorize the ancestry of the individual? If you analyze this origin of population X you should make a reference without population X and analyze all the individuals from the population with all the surrounding populations in the reference. However if you are interested in origin of a single individual, than the population X got to be in reference as well.

If your result will be that 90% of SNPs of your individual correspond to your population X, then it means that's a normal X fellow. Why it should be a problem?

-- edit --

I just found this preprint that highlights potential pitfalls in interpretation of STRUCTURE/ADMIXTURE results.


Improve a reference genome with sequencing data

One approach to this is to use whatever data you have to iteratively update the reference genome. You can keep chain files along the way so you can convert coordinates (e.g. in gff files) from the original reference to your new pseudoreference.

A simple approach might be:

Align new data to existing reference
Call variants (e.g. samtools mpileup, GATK, or whatever is best for you)
Create new reference incorporating variants from 2
Rinse and repeat (i.e. go to 1)

You can track some simple stats as you do this - e.g. the number of new variants should decrease, the number of reads mapped should increase, and the mismatch rate should decrease, with every iteration of the above loop. Once the pseudoreference stabilises, you know you can't do much more.


Do I need to remove water molecules before protein-ligand docking?

Yes, if you are using implicit solvent you ought to remove all of them.
Exception. Unless the crystallographic waters are trapped essential structural components of your active site and your ligand is not going to replace it, in which case keep it because its remove changes the energy of the active site —say you have a catalytic triad, with a water that is pinned between the acid (e.g. aspartate) and the base histidine, its removal would change the tautomerisation of latter and result in the nucleophile being protonated (and not nucleophilic). Having said that, if you don't have time to look at your structure (or cannot because of automation) and are okay with some artefacts from incorrect protonation strip them.
If you are not using implicit solvent, then it depends on the tool used —most require them stripped and re-added. If you do keep them, make sure that the crystallographic water (HOH) has the same residue code as the added water (TP3, WAT etc.).


How to split bedfile into non-overlapping regions and compute aggregation function on duplicate segments

To make disjoint intervals, you could use BEDOPS bedops --partition, piping to bedmap --mean to get the mean signal over disjoint regions.

Starting with the input bedgraph file, convert it to five-column BED with GNU awk, putting the signal in the fifth column per UCSC convention:

$ awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }' /tmp/in.bedgraph | sort-bed - > /tmp/in.bed
$ cat in.bed
chr9    8855925 8986425 .   -1
chr9    8855965 8986325 .   0
chr9    8893857 9022105 .   -1

The disjoint set from these intervals would look like this:

$ bedops --partition /tmp/in.bed
chr9    8855925 8855965
chr9    8855965 8893857
chr9    8893857 8986325
chr9    8986325 8986425
chr9    8986425 9022105

We can pipe these intervals to bedmap --mean with the original five-column BED file as a "map" file, calculating the mean signal over the disjoint intervals:

$ bedops --partition /tmp/in.bed | bedmap --echo --mean --delim '\t' - /tmp/in.bed > /tmp/answer.bedgraph
$ cat /tmp/answer.bedgraph
chr9    8855925 8855965 -1.000000
chr9    8855965 8893857 -0.500000
chr9    8893857 8986325 -0.666667
chr9    8986325 8986425 -1.000000
chr9    8986425 9022105 -1.000000

The --mean option calculates the arithmetic mean. Other aggregation functions are available to process the signal or score (fifth) column in map BED files, like --stdev, --median, --mad, etc.

When calculating statistics, the default precision is six digits. The --prec <n> option can be used to specify desired precision in the score result.

See bedmap --help or the online documentation for a full description.

To avoid making and cleaning up intermediate files, bash process substitutions could be used to generate a one-liner:

$ bedops --partition <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') | bedmap --echo --mean --delim '\t' - <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') > answer.bedgraph

This can be useful for cleaner pipelines or with processing very large, whole-genome scale inputs on systems with limited disk space.


Is it possible to convert BAM file from one genome assembly to the other?

It is possible, with tools like CrossMap, but my recommendation is to re-align. Assemblies can change a lot between versions and your reads will undoubdedly align differently. Re-aligning is cheap, it will take a few days and you're good. Re-doing your analysis because you realize the coordinate translation messed something up will cost you much more time and energy.


How to build a BAM header file with htslib in C++?

You're just looking in the wrong header file. The BCF stuff will be in vcf.h (you'll see the two functions you mentioned there).
With BAM files one doesn't typically need to add/remove headers, which is why there are not convenience functions for them. To add a header:

Increment hdr->ntargets
realloc(hdr->target_len, ...) and hdr->target_name.
You likely want to append the name/length to the end of hdr->target_name and hdr->target_len, respectively.

I don't think you need to update the dictionary or plain text for most uses (or at least I've never needed to).
Removing a chromosome would follow the reverse instructions, (you'd need to shift everything in the name and length arrays if you remove entries in the middle.
WARNING, doing this is prone to error!
I would strongly cautious you to be very careful when doing this. Please remember that the chromosome that entries are mapped to isn't actually stored in alignment entries themselves. Rather, they simply have an index into hdr->target_name, which saves space. Consequently, if you start adding/removing headers to/from a BAM file with alignments then you're going to have to modify all of the alignments to keep the correct alignment:chromosome associations. Otherwise you'll end up with completely useless results (this is the same caveat as with samtools reheader).


Receptor-Ligand R-package for visualization

Two nice tools for receptor-ligand interaction analysis are:

CellPhoneDB by the Teichmann lab, although this is not an R package...
NicheNet by the Saeys group

Note that these tools are mainly meant for analysis of receptor-ligand interactions, although I think NicheNet also has some visualization options. CellPhoneDB needs to be run from a shell, but the results can easily be imported into R for visualization. You might check the papers on these 2 methods for some ideas of visualization. Combining Seurat's visualization options and ggplot2 should get you already pretty far.


How can I read FCS files using open source libraries?

A few years ago, I wrote a python script to convert FCS files into tab-separated format. It was far from handling all the possibilities that the format description offers, but at least it worked for some of the files produced on one of our machine: http://www.igh.cnrs.fr/equip/Seitz/en_equipe-programmes.html

The format documentation I found enabled decoding (see section 3 of the pdf you mention), but it requires reading data in binary mode.

The general idea of this format (and, I guess, many other binary formats), is that there is a header zone at the beginning of the file with a defined number of fields encoding numbers indicating how the rest of the file is structured. So a first phase is to parse this header, following the description given in the documentation of the format. The information extracted from the header tells where to find the data and how it is encoded, still according to rules described in the format documentation.

In case this may be useful, and for the record, here is the code from the above-mentioned script (after stripping comments, some of which are merely copied from the format documentation, and adding a few ones):

#!/usr/bin/env python
"""This script tries to read FCS flow cytometry data.
Format parsing inspired by information found here:
http://isac-net.org/Resources-for-Cytometrists/Data-Standards/Data-File-Standards/Flow-Cytometry-Data-File-Format-Standards.aspx
"""

import re
# To decode binary-encoded data
import struct
import sys

class Parameter(object):
    """This object represents one of the parameter types that are present in a DATA segment of a FCS file."""
    __slots__ = ("p_name", "p_bits", "p_range", "p_ampl", "parser")
    def __init__(self, p_name, p_bits, p_range, p_ampl):
        self.p_name = p_name
        self.p_bits = p_bits
        self.p_range = p_range
        self.p_ampl = p_ampl
        # Function for parsing a value of the parameter in the data segment
        self.parser = None

##############################################
# Here starts the parsing of the header part #
# which tells where the other parts are.     #
##############################################

f = open(sys.argv[1], "rb")
# The format name is encoded in 6 letters
# An ASCII letter is coded with one octet
file_format = "".join([f.read(1) for __ in range(6)])
sys.stdout.write("Format: %s\n" % file_format)
# The format descriptions reserves 4 octets that we skip
skip = f.read(4)
# 8 octet chunks encode the start and end positions
# of different parts of the data
text_start = int(f.read(8).strip(" "))
text_end = int(f.read(8).strip(" "))
data_start = int(f.read(8).strip(" "))
data_end = int(f.read(8).strip(" "))
analysis_start = int(f.read(8).strip(" "))
analysis_end = int(f.read(8).strip(" "))

if (analysis_start and analysis_end):
    sys.stderr.write("Cannot deal with ANALYSIS segment of an FCS file.\n")

####################################################
# Here starts the parsing of the "TEXT" portion    #
# which describes how the data proper is organized #
####################################################
f.seek(text_start)
# The first character in the primary TEXT segment is the ASCII delimiter character.
sep = f.read(1)
if sep not in ["_", "@"]:
    alt_sep = "_@_"
elif sep not in ["_", "|"]:
    alt_sep = "_|_"
else:
    assert sep not in ["+", "|"]
    alt_sep = "+|+"
text_segment = f.read(text_end - text_start)

fields = text_segment.split(sep)

info = {}

i = 0
while i < len(fields) - 1:
    key = fields[i]
    i += 1
    val = fields[i]
    i += 1
    # Keywords are case insensitive, they may be written in a file in lower case, upper case, or a
    # mixture of the two. However, an FCS file reader must ignore keyword case. A keyword value may
    # be in lower case, upper case or a mixture of the two. Keyword values are case sensitive.
    info[key.upper()] = val
print "%s events were detected." % info["$TOT"]
print "Each event is characterized by %s parameters" % info["$PAR"]

if info["$NEXTDATA"] != "0":
    sys.stderr.write("Some other data exist in the file but hasn't been parsed.\n")

# L - List mode. For each event, the value of each parameter is stored in the order in which the
# parameters are described. The number of bits reserved for parameter 1 is described using the
# $P1B keyword. There can be only one set of list mode data per data set. The $DATATYPE
# keyword describes the data format. This is the most versatile mode for the storage of flow
# cytometry data because mode C and mode U data can be created from mode L data.
assert info["$MODE"] == "L"

parameters = []

# indices of the parameters
p_indices = range(1, int(info["$PAR"]) + 1)
for i in p_indices:
    p_name = info["$P%dN" % i]
    p_bits =  info["$P%dB" % i]
    p_range = info["$P%dR" % i]
    p_ampl =  info["$P%dE" % i]
    parameters.append(Parameter(p_name, p_bits, p_range, p_ampl))

sys.stdout.write("The parameters are:\n%s\n" % "\t".join([par.p_name for par in parameters]))

# How are 32 bit words organized
if info["$BYTEORD"] == "4,3,2,1":
    endianness = ">"
else:
    endianness = "<"
    assert info["$BYTEORD"] == "1,2,3,4"

    # I stripped a long comment which is just a copy of the documentation
# Type of data:
if info["$DATATYPE"] == "I":
    for par in parameters:
        nb_bits = int(par.p_bits)
        assert nb_bits % 8 == 0
        nb_bytes = nb_bits / 8
        # Determine format string for unpacking (see https://docs.python.org/2/library/struct.html)
        if nb_bytes == 1:
            c_type = "B" # unsigned char
        elif nb_bytes == 2:
            c_type = "H" # unsigned short
        elif nb_bytes == 4:
            c_type = "L" # unsigned long
        elif nb_bytes == 8:
            c_type = "Q" # unsigned long long
        else:
            raise ValueError, "Number of bytes (%d) not valid for an integer (see https://docs.python.org/2/library/struct.html#byte-order-size-and-alignment)." % nb_bytes
        fmt = "%s%s" % (endianness, c_type)
        p_range = int(par.p_range)
        def parser(data):
            value = struct.unpack(fmt, data.read(nb_bytes))[0]
            try:
                assert value < p_range
            except AssertionError:
                print "Value %s higher than %d" % (str(value), p_range)
            return value
        par.parser = parser
    pass
else:
    raise NotImplementedError, "Only the parsing of integer value has been implemented so far."

out_file = open(sys.argv[2], "w")
out_file.write("#amplification_types\t" + "\t".join([par.p_ampl for par in parameters]) + "\n")
out_file.write("parameters\t" + "\t".join([par.p_name for par in parameters]) + "\n")
i = 1
##############################################
# Here starts the parsing of the data proper #
##############################################
f.seek(data_start)
while f.tell() < data_end:
    values = []
    for par in parameters:
        values.append(par.parser(f))
    out_file.write("%d\t" % i + "\t".join(map(str, values)) + "\n")
    i += 1
out_file.close()
f.close()


Calling for one specific SNP from multiple sequencing runs

Getting SNPs from whole-genome data is usually done via mapping sequencing reads on the reference (bwa-mem or bowtie are two popular mappers) and then using one of variant calling tools (like freebayes, samtools, or GATK). However, I would also look up if there is something virus-specific.
According to this answer, you could get variant calls of a specific region using samtools if you give it a bed file that contains the region of interest (three columns chr    start    stop format). For example using human reference
samtools mpileup -uDl regions.bed -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - > RAL_samtools.raw.bcf

You could also specify a bunch of files, not just two and it would do the variant calling on them all.


Question 1: Is there a tool for imputing phenotypes while also using genetic data?

Yes, try Andy Dahl's Phenix from the Marchini lab.

Disclosure: I did my PhD in the Marchini lab, although I was not directly involved in this research.

Question 2: Is there a way to impute categorical data using methods intended for continuous data?

If the variable of interest has only two categories, then you should be able to simply recode the categories as arbitrary (differing) numbers and get sensible results. But do make sure you understand the tool you are using to see if this will work.

For more than two categories, if your imputation approach allows the simultaneous imputation of multiple variables, you could try recoding your categories as contrasts or one-hot.  After recoding, you could impute the variables resulting from recoding together.


I would strongly discourage you from making discontinuous axis, it's going to be very confusing for a reader.

The facet plot you proposed seems like a good solution to me. Alternatively you can use log transformation. To demonstrate I made it on simulated data that look appox like yours :

set.seed(940401)

data.plot <- data.frame(Methyl_Average = c(0.05, rbeta(699, 2, 5) + 0.15, runif(300, 0.18, 0.85)),
                        Average_Gene_Expression = c(165, exp(rnorm(999, 3))))

library(ggplot2)
require(gridExtra)

g <- ggplot(data.plot, aes(x = Methyl_Average, y = Average_Gene_Expression))
g1 <-  g + geom_point(color = "steelblue") + theme_bw(base_family = "Times") +
       ggtitle( "Scatter Plot") +
       theme_bw(base_family = "Times") +
       theme(plot.title = element_text(hjust = 0.5))

sc_plot <- g1 +
           labs(x = "Differential Methylation", y = "Differential Expression")
logsc_plot <- g1 +
              labs(x = "Differential Methylation", y = "Differential Expression (log10 scale)") +
              scale_y_continuous(trans='log10')

grid.arrange(sc_plot, logsc_plot, ncol = 2)


Parallel version of cn.mops' bamDataRanges

This is an option for the function documented in its help page. Just run ?bamDataRanges and you will see:

Usage:

     getReadCountsFromBAM(BAMFiles, sampleNames, refSeqNames, WL = 25000,
       parallel = 0, ...)

[...]

parallel: The number of parallel processes to be used for this
          function. Default=0.

So all you need to run with 10 cores is:

bamDataRanges <- getReadCountsFromBAM(BAMFiles, parallel=10)


What does genotyping data and beta values for SNPs mean for GWAS study?

The beta statistic is a common output field produced by programs that carry out genome-wide associations. There's a good explanation here.
I interpret a beta statistic as the degree to which a change in genotype from the reference allele alters the phenotype of interest. The association (and beta statistic) can be negative or positive. It should always be tied to a particular variant, usually called the "risk allele".
The beta statistic is related to the number of individuals found in each group with a particular genotype. In that sense it could be argued to be related to the population size (at least, assuming pre-grouping of populations by genotype), but that's not a typical interpretation.
The GWAS catalogue is a good place to start off looking for data associated with GWAS:
https://www.ebi.ac.uk/gwas/home
However, I think that's for summary statistics. If someone has asked you for raw genotype data, that is trickier to get, and tends to depend on the condition. Genotype data is personally-identifying information, and there are lots of issues associated with a blanket public release of such information, which is why it's so hard to find.
Without more information on the particular phenotype, I can't give much guidance, but useful places to look for most traits are dbGaP and UK Biobank. Both of those data sources need forms filled out and approved (to demonstrate an ethically-responsible research aim) in order to access genotype-level data.
There are a few public datasets around that include data uploaded by individuals (who have explicitly consented to their genotype data being public), but it tends to be trait-specific.


How to get results from Homo.sapiens package in bioconductor for a specific reference

I couldn't find how is this package build, my guess is that is as in the example of the OrganismDBI package. See that vignette section to build your own data package with the genome you want.

The locations are different because you ask for the transcript start and end, if you add the "TXNAME", you will (probably) see there are different transcripts, so different mappings are reported (untested).

Update

The example of the OrganismDBI package builds correctly the Homo.sapiens package, which appears to be linked with hg19. I provide here the relevant code to build the package for hg19 and hg38.

> source("https://bioconductor.org/biocLite.R")
> biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
> biocLite("TxDb.Hsapiens.UCSC.hg38.knownGene")

> library(OrganismDbi)

# instructions for hg19, for hg38 change accordingly
> gd <- list(join1 = c(GO.db="GOID", org.Hs.eg.db="GO"),
             join2 = c(org.Hs.eg.db = "ENTREZID",
                       TxDb.Hsapiens.UCSC.hg19.knownGene = "GENEID"))
> destination <- tempfile()
> dir.create(destination)
> makeOrganismPackage(pkgname = "Homo.sapiens.hg19", graphData = gd,
                      organism = "Homo sapiens", version = "1.0.0",
                      maintainer = "Maintainer<maintainer@email>",
                      author = "Author Name", destDir = destination,
                      license = "Artistic-2.0")

Creating package in /var/folders/sf/tfrswpc96jl8b4520mn784440000gn/T//RtmpzR4YnY/file111b490af70f/Homo.sapiens.hg19

# installing package
> install.packages(destination, repos = NULL, type="source")

Then, you can load and use the packages with:

> library(Homo.sapiens.hg19)
> library(Homo.sapiens.hg38)

> select(Homo.sapiens.hg19, keys=c('KRAS'),
     cols = c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'),
     keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863
> select(Homo.sapiens.hg38, keys = c('KRAS'), cols =
       c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25204789 25250931
2   KRAS   chr12 25209168 25250936
3   KRAS   chr12 25209431 25250803
4   KRAS   chr12 25233819 25250929


Does UCSC definition for telomere making sense?

It seems so, yes.

UCSC uses a 10Kb definition and representation for telomeres in human. In the table you shared this is consistent across all chromosomes, not only chr21. This 10kb representation is an approximation to the mean value observed in newborns and is supported by evidence like the one linked below.

Taken from Wikipedia:

A telomere is a region of repetitive nucleotide sequences at each end of a chromosome, which protects the end of the chromosome from deterioration or from fusion with neighboring chromosomes. For vertebrates, the sequence of nucleotides in telomeres is AGGGTT, with the complementary DNA strand being TCCCAA, with a single-stranded TTAGGG overhang. This sequence of TTAGGG is repeated approximately 2,500 times in humans. In humans, average telomere length declines from about 11 kilobases at birth[ to less than 4 kilobases in old age, with the average rate of decline being greater in men than in women.

2,500 x 6 (TCCCAA) = 15,000
(approx. 10K)

Two interesting additional facts:

The minimum telomere length needed to ensure human telomere protective stability in white blood cells is 3.81 kb
Mouse has a 100Kb telomere representation in USCS


Although I've never tried it myself, I think what you need can be achieved through the LogCombiner program from the BEAST package -- assuming the species trees files are comparable (ran with same parameters, same species names, etc.) I am also assuming that you meant 'five species tree files', each file with the posterior sample of species trees.

However, I do believe you should try to run all loci within the same analysis, or at least do some prior check to see if the distributions between data sets are comparable. LogCombiner assumes independent runs from the same data.

I would also recommend taking a look (and following) the BEAST google group (which includes starBEAST discussions).


How drastic would the result be if I use fasta or reference assembly from ucsc and gtf from gencode?

Never use genomes or annotations from UCSC, they're poorly versioned and only recently with mouse and human have they even included all of the contigs. For fasta/GTF files from early in the GRCh38 release, you can tell whether you're using UCSC or Gencode by the presence/absence of _random contigs, which will only exist for UCSC. These were mostly later split into the actual contigs, so recent download from UCSC should more closely match what you find at Gencode/Ensembl. Further, that time predated UCSC beginning to adopt Gencode's vastly superior annotations, so if your GTF file has instances where the same gene ID is on either multiple strands or multiple chromosomes (this is obviously biologically impossible) then you have a UCSC GTF file.

In general, with early GRCh38 releases your only real issues will be with _random contigs, which are a minority of the genome and don't have all that many genes. But really, you should be keeping track of the sources of your files and ensuring that they're compatible.

Update: I should expand a bit on my "Never use genomes or annotations from UCSC" comment. In point of fact the genomes themselves aren't so terrible. Early on UCSC had the bad habit of concatenating contigs together into _random "chromosomes", but they seem to have mostly kicked that habit as of late. Note, however, that there are no versions for their genomes. Since reference genomes continue to get updates over time (mostly through the addition of patches) the lack of actual versions means you have to manually check if a recently downloaded file matches what may have been downloaded either previously or by someone else. This has obvious consequences when it comes to reproducibility. The same issue occurs for annotations from UCSC, but they have the additional problem of historically having biologically incoherent concepts of genes. That is, they will contain the same gene in multiple places with multiple orientations, which will break many tools in both obvious and completely unclear ways. For example, DEXSeq will simply break with an error message if given a UCSC annotation, since they break biological plausibility. If you were to use these annotations files with deepTools, you wouldn't get an error message, but the resulting output would be only partial, due to the biologically impossible annotation effectively corrupting most obvious ways of storing annotation data in a data structure (i.e., you can no longer treat IDs as unique). This could have downstream ramifications on biological interpretation of results.


Where to Download Cancer Raw Reads (fastq)?

You can try to browse SRA (sequence read archive) at https://www.ncbi.nlm.nih.gov/sra.
To download fastq files from SRA, you'll need their sratoolkit: https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software.

The easy command to download fastq files by their accession number is :

fastq-dump --split-files accession_number


The only reason to normalize for GC content in RNA-seq is if it differs notably between samples/groups. If that's not the case and you aren't trying to compare genes withing samples then you have no reason to try to account for GC content. The same goes for ATAC-seq, though there the danger with trying to normalize for GC content is that you then mask changes in global accessibility.
In short, normalize for things that should cause problems with your analysis and nothing else.


There are three possible problems that come to mind.

Masking
Blast will mask low complexity regions by default. Since your sequence is nothing but Gs, it is a safe bet that it is being masked, so no hits will be found for it.

Score/e-value thresholds
Another source of complication is that even if a match is found, that match will have very bad scores. Both the actual score of the alignment and the e-value will be very bad. Since this is such a simple sequence, it will always score badly.

Word-size
The way blast works is (simplifying a bit) by finding a match for N residues (the word size) and then extending that match if extending increases the score. If your query sequence is shorter than the word size, no match will ever be found.

With all that said, the specific problem in your case seems to be a combination of the word-size and masking. If you run your query setting the word-size to 11 and disabling masking, as shown below, you should get results.
blastn -db A.fna -query B.fna -word_size 11 -dust no

That should return two hits (two because the word size is one less than the query size, so it can align perfectly at two locations):
> A
Length=64

 Score = 23.3 bits (12),  Expect = 5e-05
 Identities = 12/12 (100%), Gaps = 0/12 (0%)
 Strand=Plus/Plus

Query  1   GGGGGGGGGGGG  12
           ||||||||||||
Sbjct  11  GGGGGGGGGGGG  22

 Score = 23.3 bits (12),  Expect = 5e-05
 Identities = 12/12 (100%), Gaps = 0/12 (0%)
 Strand=Plus/Plus

Query  1   GGGGGGGGGGGG  12
           ||||||||||||
Sbjct  12  GGGGGGGGGGGG  23

Note that the default word size for BLASTn searches is usually 11. I don't know why the blast+ version (2.2.31+) I have on my Ubuntu machine didn't use that as a default. The manual says it should. But I got no hits when I didn't explicitly set the word size to 11.


How to modify DNA evolution model to fit actual data?

Mutation is not a random process because there are many mutational biases, e.g. the transition vs. transversion ratio. Ideally use a proper evolutionary model. Maybe try a tool that incorporates substitution and indel models. I'm not up-to-date on the best packages, but e.g.

INDELible http://abacus.gene.ucl.ac.uk/software/indelible/
Evolver https://www.drive5.com/evolver/ (not to be confused with the PAML package of the same name)


I'll fix the bamCoverage issue in the 2.6.0 release (it may already be fixed since I rewrite how bigWig files are made for that release). You can follow my status on the issue page. In principle you should be able to do this, it's just a deepTools limitation.

This should work fine with computeMatrix, though keep in mind that you're just going to get a big block of NA values.

Update: This error ends up coming from pyBigWig. I've just found that bedGraphToBigWig also can't make a bigWig file with now entries, so I suspect this is sort of undefined behavior. Regardless, I'll alter pyBigWig to produce a file that it can open and use in a meaningful manner.


Produce .vcf file of ALL mutations in .bam file

Ensure that the variant caller stringency is relaxed all the way. For example in freebayes use --min-alternate-count 1 --min-alternate-fraction 0 .


Where can I find Single Cell Data with Location "Coordinates"?

Does single cell data typically have the following meta-data: the
"coordinates" (e.g. on a tissue, adjacent tissues) saying where each
cell in the sample was located relative to other cells?

No. Typical scRNA-seq is just capturing random cells in a tube with no additional information. The technology you are looking for is spatial transcriptomics where you are measuring RNA levels at a particular location, but that is still not on the single-cell level.

If not, is it possible to reconstruct this with other meta-data on the
cells?

There have been efforts to reconstruct spatial relationship. For example: Satija et al. 2015. There should be more recent approaches that I am not aware of.


Presuming you still have your BAM file sitting around somewhere:

#!/usr/bin/env python
import pysam
import pyBigWig

bam = pysam.AlignmentFile("alignments.bam")
ifile = pyBigWig.open("your current file.bw")
ofile = pyBigWig.open("your new file.bw", "w")

# Create/Add the header
hdr = []
for chrom in bam.header['SQ']:
    hdr.append((chrom['SN'], chrom['LN']))
bam.close()
ofile.addHeader(hdr)

# Write the bigWig
for _ in hdr:
    try:
        ints = ifile.intervals(_[0])
        starts = [x[0] for x in ints]
        ends = [x[1] for x in ints]
        vals = [x[2] for x in ints]
        ofile.addEntries([_[0]] * len(starts), starts, ends=ends, values=vals)
    except:
        # intervals() will throw an error if the chromosome isn't present
        pass
ofile.close()
ifile.close()

The reason you can't append to a bigWig file is that it has a various indices and zoom levels typically built in. You'd have to add to them without overwriting neighboring blocks of data. In theory this might sometimes be theoretically possible, but not typically and I'm not aware of tools that support this.


Seurat CCA, AlignSubspace, Argument dist.method is only useful with multivariate timeseries

I don't have enough reputation to comment so I will briefly answer here

This warning is fixed in Seurat 2.3.2. Without upgrading you can ignore this message. It doesn't mean anything for the actual analysis performed so your results won't be flawed.


How to get rid of the temp files while using "featureCounts" for extracting readcounts from bam files?

You shouldn't delete them while featureCounts is running, it's creating read-name-sorted temporary files so it can more easily look at both mates together. Just wait for it to finish running and it will clean up the files itself. Doing so beforehand is liable to cause an error.


What is the best way to programmatically convert Ensembl ids from one release to other? e.g. transcript ids from release 75 to 100

The ID history converter is available as a Perl script.
This script accesses these MySQL tables, so you could query these tables directly.
Note that the mapping session is only one release at a time, so you would need to run 25 queries to go from release 75 to 100.


Is there an efficient way to check an input BAM in R?

You won’t know whether it’s a valid BAM file until you attempt reading it in full; even Devon’s method (which I recommend) only does a superficial check — i.e. it checks whether (parts of) the BAM header can be read correctly. However, this doesn’t guarantee that the file isn’t truncated or contains invalid data.

This is of course a more general problem, but it’s perhaps a bit stranger for BAMs (and other biological formats), as it’s relatively frequent to end up with technically invalid files that were produced by existing software. The only recourse is actually trying to use the file, and to handle errors appropriately.

For a very quick check™ you could also pass the file through gunzip and test whether the “magic bytes” at the beginning of the file equal 'BAM\1':

bgzf_stream = gzfile(filename, 'r')
magic = readChar(bgzf_stream, 4)
stopifnot(identical(magic, 'BAM\1'))


Annotating gene names or gene IDs to a dataframe containing SNPs?

I assume your SNP file contains genomic positions - please clarify if not. For example, you could have an rsID column, chromsome, and position column that contain "rs3088379", "chr1", "54225281".
Approach #1: use R, as you suggest.
ensembldb is an R package that may be what you're looking for: https://bioconductor.org/packages/release/bioc/html/ensembldb.html.
You could map the genomic positions of your SNPs to gene names by adapting their workflow: https://bioconductor.org/packages/release/bioc/vignettes/ensembldb/inst/doc/coordinate-mapping.html.
After you've mapped the genomic position, you could filter for your gene name.
Approach #2: use bedtools.
You can convert your Excel file of genomic positions to a tab-delimited BED format. https://bedtools.readthedocs.io/en/latest/content/overview.html has an overview. So, your first SNP would turn into
chr1 54225280 54225281 rs3088379

--> make sure you turn 1-based to 0-based coordinates!
You could then download the coordinates of your gene in BED format from your favorite source, like UCSC Table Browser.
Then, you run bedtools intersect with your SNP BED and your gene BED. What remains is just that gene's SNPs.


Gene Ranking - signal to noise ratio used in GSEA-P algorithm?

The noise is comming from the subset.mask, which is created above in a loop with the number of permutations.

for (r in 1:nperm) { #L90
....
        subset.mask[, r] <- as.numeric(c(subset.class1, subset.class2)) # L107

So by multiplying the random selection of subsets by the expression we get the "noise". Later on line 251 we get the ratio of noise/signal by dividing the signal from the contrast and the signal by the noise.

Those transformation on the middle seem to come from normalizing the signal in order to be able to compare the resampling, but I don't fully understand why they do this. Also note that the R scripts are outdated and no longer maintained, so they can contain errors.


I've given an answer about using how I normalise for PCA here. PCA works best if the data are normal-ish, but as models go it is fairly robust to fairly large deviations from normal. Transforming is not essential, but if data can be easily transformed into something that looks a bit more normal, then conclusions derived from looking at the PCA should be a bit more robust.

I'm not familiar with canonical-correlation analysis (CCA), but because I see sums and square roots in the formulas, I would assume it has similar data cleaning requirements to PCA. Unfortunately the R documentation for cancor doesn't have any information about any assumptions made on the data.


Ensembl to PDB

So the dataset you need is from SIFT. There is one with Ensembl IDs. Majorly this also includes what chain is your peptide and what is the offset. About the latter, I actually wrote a blog post about it recently.

PDB download

Using Python, the Biopython PDB module is very good. Although the PyMol module is more powerful. But as you want only coordinates (right?), you'd best use the former, which has a better documentation.


Why am I getting empty expression data from GEO?

I am sure the data is legit, you are just approaching it incorrectly. getGEO is an application for microarray data, not for digital count data such as (sc)RNA-seq, therefore what you aim to do is simply not possible by design. Unless you want to start from the raw reads why not taking the at the file named GSE116237_scRNAseq_expressionMatrix.txt.gz provided at the bottom of https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116237?


How to measure the total size of a fastq file in base pairs?

I've been using this:
cat file.fastq | paste - - - - | cut -f 2 | tr -d '\n' | wc -c
Explanation :

paste - - - - : print four consecutive lines in one row (tab delimited), to merge the info for each read

cut -f2 : print only the second column, to access the sequence after the paste

wc -c : count the characters

tr -d '\n': to remove from count the eventual newline characters

(a tip for your googling: try search for "counting number of bases in fastq file")


GRCm37-designed exon target enrichment, which reference to use?

The only possible issue with using GRCm38 as the reference is that since GRCm37 was used in the bait design there are likely a few probes that will show artificial structural variations, due to errors in GRCm37. However, these will be present in all of your mice and if you use liftOver to look at the regions around the probes then such regions should become obvious. In fact, that'd be a convenient thing to do if you want to keep your (likely better) GRCm38 alignments, flank each probe site by a few hundred bases and then use liftOver to get the GRCm37 coordinates. If you go through the result you'll be able to find the false-positive structural variants due to the baits being designed with an older reference.

Alternatively, just use GRCm37 and hope that any issues with it don't result in false-negative SV calls.


Is it possible to perform MinION sequencing offline?

Short answer: yes, but you need to get permission (and modified software) from ONT before doing that.

... but that doesn't tell the whole story. This question has the potential to be very confusing, and that's through no fault of the questioner. The issue is that for the MinION, sequencing (or more specifically, generating the raw data in the form of an electrical signal trace) is distinct and separable from base calling. Many other sequencers also have distinct raw data and base-calling phases, but they're not democratised to the degree they are on the MinION.

The "sequencing" part of MinION sequencing is carried out by ONT software, namely MinKNOW. As explained to me during PoreCampAU 2017, when the MinION is initially plugged into a computer it is missing the firmware necessary to carry out the sequencing. The most recent version of this firmware is usually downloaded at the start of a sequencing run by sending a request to ONT servers. In the usual case, you can't do sequencing without being able to access those servers, and you can't do sequencing without ONT knowing about it. However, ONT acknowledge that there are people out there who won't have Internet access when sequencing (e.g. sequencing Ebola in Africa, or metagenomic sequencing in the middle of the ocean), and an email to <support@nanoporetech.com> with reasons is likely to result in a quick software fix to the local sequencing problem.

Once the raw signals are acquired, the "base-calling" part of MinION sequencing can be done anywhere. The ONT-maintained basecaller is Albacore, and this will get the first model updates whenever the sequencing technology is changed (which happens a lot). Albacore is a local basecaller which can be obtained from ONT by browsing through their community pages (available to anyone who has a MinION); ONT switched to only allowing people to do basecalling locally in about April 2017, after establishing that using AWS servers was just too expensive. Albacore is open source and free-as-in-beer, but has a restrictive licensing agreement which limits the distribution (and modification) of the program. However, Albacore is not the only available basecaller. ONT provide a FOSS basecaller called nanonet. It's a little bit behind Albacore on technology, but ONT have said that all useful Albacore changes will eventually propagate through to nanonet. There is another non-ONT basecaller that I'm aware of which uses a neural network for basecalling: deepnano. Other basecallers exist, each varying distances away technology-wise, and I expect that more will appear in the future as the technology stabilises and more change-resistant computer scientists get in on the act.

Edit: ONT has just pulled back the curtain on their basecalling software; all the repositories that I've looked at so far (except for the Cliveome) have been released under the Mozilla Public License (free and open source, with some conditions and limitations). Included in that software repository is Scrappie, which is their testing / bleeding-edge version of Albacore.


What you are looking for is SNP annotation. If you have the chromosome:position reference and alternate alleles for your SNPs of interest, it can be as simple as uploading them to the variant effect predictor.

http://grch37.ensembl.org/Homo_sapiens/Tools/VEP

This will give you the predicted protein change and novelty of the variant with respect to known variation in large scale sequencing projects like the 1000 genomes.

If you cannot afford to upload data on the web, there are downloadable command line software like SNPeff, Variant effect predictor etc but it would require some more bioinformatics knowledge of converting your data to a VCF format and obtaining a human reference genome.


What does an FDR value of 1 in RNA-seq mean?

FDR stands for False Discovery Rate. It is a statistic tool used in multiple hypothesis testing.

As you may know, when you use a p-value cutoff (usually 0.05) for your experiments, it means strictly speaking that "is there was actually no signal, there would be a probability of 0.05 to observe this kind of extreme values". This can be understood as "there's a 5% chance that what you call a hit is actually not signal, just some extreme value that happened by chance, and really there's nothing to see here".

This definition has deep consequences if you run many experiments: imagine that you run 100 experiments (say, you test the effect of 100 different molecules), if none of them has any effect, on average you could still get 5% of false positives. So if you run these 100 experiments, and you get, say, 10 hits with a p-value below your thresold of 0.05; how can you tell which ones are real and which ones are just chance in the absence of biological signal?

This is where you can use FDR.

If you control for a false discovery rate of, say, 0.2 (20%), it means that after a few computations (the Benjamini-Hochberg procedure described on wikipedia is simple and commonly used) you will lower your p-value cutoff so as to be sure that, among the N hits (probably less than the 10 you initially selected above) that you got above, a maximum of 20% are actually false positives (i.e. no biological effect, but surprisingly high/low signal values by pure chance).

In the article you mention they have a different FDR value for each sample, which is unusual. The values are not round numbers, so it seems likely that they are numbers obtained after applying the FDR correction (i.e. they are "transformed p-values") on which they based their decision to include the samples or not. The FDR cutoff that they chose should be written somewhere in the main text.


How to calculate overlapping genes between two genome annotation versions

Short Answer:
In my opinion, my approach would be to pull out the CDS exons and run bedtools on those.

A Few More Details:
When you pull out the exons, make sure that you assign them all IDs if the don't already have them assigned and record which IDs "belong" to which genes. Now when you get exons that overlap, you know that they are coding and you can tie them back to which genes they originate from.


Why do molecular generation models maximize “penalized logP” as a measure of drug-likeliness?

A positive logP is good (hydrophobic), but high results in ADMET problems (cytochromes break overly hydrophobic compounds) and aggregation (gunk, tar, crud are words happily thrown around). However, in most cases the logP is not greater than +3 because good binding is result of hydrogen bonds etc —some bad algorithms will suggest a polyaromatic hydrocarbon that makes lots of pi-pi interactions but these are few.
Confusingly, a high synthetic accessibility score is bad (hard). 3 is the average and 10 is basically something impossible. This score is how hard it is to make (i.e. pricier), which is one of the major drivers.
The number of cycles is number of rings. Even if you get pi-pi, cation-pi and sulfur-pi bonds you also get aggregation due to stacking (a crud). So is somewhat bad.
So basically that score is
high=good minus high=bad minus high=bad

If difference in relative Gibbs energy of binding is thrown it, negative ∆∆G is good, so that would be subtracted. As its atom count dependent ligand efficiency is better, but that ranges from 0 to -1.2 kcal/mol. So would need scaling.


Something like this book chapter might help on understanding FDR. The packages or functions you call differ in their estimation of pi0 or the proportion of null hypothesis.
Basically, the p.adjust(..method = "BH") is a more conservative method with Benjamini-Hochberg, you are assuming the proportion of null features to be 1.
When you use fdrtool, you are estimating the proportion of null hypothesis and from there, the local FDR. This in theory gives you more power.
However, I have to note that in your code, you are not using fdrtool to estimate the fdr as explained above, because:
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = FALSE)

With this, the raw p-values are stored in FDR.ddsRes$pval and the q value, which is the corrected p-value, is stored in FDR.ddsRes$qval. In your code you are basically taking the raw p-values, and applying Benjamin-Hochberg:
p.adjust(FDR.ddsRes$pval, method = "BH")

You should apply a two tailed test. From the wald stats, it should be:
derived_p = 2*pnorm(-abs(ddsRes$stat))

In your code, you are running it wrongly by passing the test statistic into fdrtools, we can compare the raw p-values:
head(data.frame(stat = ddsRes$stat,derived = derived_p,
DESeq2=ddsRes$pvalue,fdrtools=FDR.ddsRes$pval))

        stat    derived     DESeq2  fdrtools
1  1.2405796 0.21476110 0.21476110 0.1216358
2  1.2314421 0.21815755 0.21815755 0.1244054
3 -0.3292516 0.74196552 0.74196552 0.6811997
4 -0.4964574 0.61957174 0.61957174 0.5356140
5 -0.4101163 0.68172066 0.68172066 0.6088427
6 -1.8007943 0.07173531 0.07173531 0.0246428

To summarize, use the adjusted p-value returned by DESeq2. The code you have above is wrong unless there's a good reason to do a one-tail test.


You just need to strip the / from $d, for which there are a number of options:

for d in */ ; do
    dname=`basename $d`
    findPeaks $d -style factor -o ${dname}.txt
done

or

for d in */ ; do
    dname=${d%%/}
    findPeaks $d -style factor -o ${dname}.txt
done

There are likely other ways one could go about this. Having said that, since you're not running these in parallel you could simply move the default output file (peak.txt) to a file name of your choice (using one of the strategies above).


What software can I use to develop a population haplotype network?

The software you are looking for is Network 5.0.0.3

Phylogenetic Network Software . Network © Copyright Fluxus Technology Ltd 1999–2018. Retrieved from https://www.fluxus-engineering.com/sharenet.htm

Whilst its old code it is still regularly used in pop gen for example there guys here recently used it. It will do the job, although personally I think Google Fusion tables are very cool, but that is complicated to use properly. R has some cool stuff in this context BTW, but in context I'd keep it simple (i.e. Network 5)


A good tool for gene locus visualization

I think in R, Gviz is the best choice. I have used it for human sequences but it will probably also be suitable for Prokaryotes.

Another option is MS powerpoint, if you only need the arrows...


Snakemake wildcard constraints reject too much

Does anyone know if there is anything special about the snakemake regexes? Should they not be normal python regexes?

No, nothing special.  The Snakemake documentation just links to the Python re module documentation.
I think this is something strange in how the anchors ^ and $ behave when matching wildcards, and not anything else in your particular regular expression.  Check out the pattern-matching with this dummy Snakefile:
rule all:
    input: ["data/subdir/example", "data/example", "example"]

rule anchored_example_of_subdir:
    """Can't add either ^ or $ or this will fail"""
    output: "data/{subdir,subdir}/example"

rule anchored_example_in_subdir:
    """Can't add ^ or this wil fail"""
    output: "data/{file,example$}"

rule anchored_example:
    """This works with both ^ and $"""
output: "{file,^example$}"

(There are no inputs or actions so it's just a matter of what filenames match what rules, and the regular expressions do nothing but match exact strings.)  This "works" as is (snakemake -nrp shows what you expect) but if you add in anchors as described in the docstrings above the upper two rules no longer match the expected filenames.  It's as if (and maybe is?) under the hood Snakemake is trying to apply the expression to the entire path.
Can you just drop your ^ and $?  The snakemake documentation implies that it's using the equivalent of re.match rather than re.search already, with wording like "Here we could restrict the wildcard dataset to consist of digits only using \d+ as the corresponding regular expression" so I don't think it's necessary anyway.

Edit: I think I found the issue in the snakemake source.  It munges your constraints together into one larger regular expression for the entire file path, so ^ and $ get embedded inside the resulting expression.  This looks like something that should at least be clarified in the documentation.


Can I use GSEA to look for enrichment in data that is grouped into sets but are not necessarily genes?

Yes, if you have a ranked list, you can use GSEA to calculate the enrichment of the strains along the ranked list.

Note that hypergeometric test is for count data, so you would test if there are more strains of one type that fit on certain conditions.

I would recommend this article comparing different enrichment methods to understand well what they do.


BUSCO does not find buscos in reference sequence of C. elegans (technical problem)

Most likely your Augustus installation is not functioning properly. Did you install Augustus from source or using conda? I believe the conda version of Augustus is known to have some issues.

Something worth trying: Running a complete busco install in a VM. This image has all dependencies installed and properly configured.


Effect of mutation in DNA sequence on transcription factor binding sites

I haven't worked on TFBSs but I really, really doubt you will be able to do this. In mutations, like in real estate, it's location, location, location. It isn't a question of how many mutations but what mutations and where. TFBSs are usually described using position weight matrices, precisely because that's how you can show the importance of different changes in different positions.

Here are a couple of random Google image search results for "transcription factor binding sites logo":

(This was taken from here and you might want to read through it.)

(source: Front. Genet., 23 February 2016 | https://doi.org/10.3389/fgene.2016.00024)

The main take home messages from these images are:

Different positions have different importance and potential to disrupt binding. In the first image, I would expect that changing the T at position -3 to anything else would be pretty bad. However, changing an A to a T at -8 would likely be neutral.
Different TFs have very different binding sites and very different binding site lengths.  Look at panel B of the second image. That one appears to be a very conservative site and likely any change would be deleterious. Conversely, changing any of the nucleotides of the 17bp spacer shown in panel A would likely not be important as long as you don't change the length.

Please note that I have not done more than glance at the sources I took these from, so my inferences may well be wrong, I only use them to illustrate the general point: you cannot really put a number to this. You will always need to look at each specific change and each specific position. In some cases, a single change can disrupt binding while in others multiple changes might have no effect at all.

You might want to look into TFBS databases like JASPAR or TRANSFAC.


Filtering bam file with bedtools intersect using strand specificity

Your bed file switch the score and strand columns. Strand should be column 6, score column 5. For example take this bed file with the 5th and 6th columns switched - zero results for strand specific intersection b/c strand is actually '.'

$ cat wrong.bed
chr1    30000000        31000000        noname  +       .
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 \
    | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 -s \
   | samtools view | wc -l
0

After fixing the bed file:

$ cat right.bed
chr1    30000000        31000000        noname  .       +
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 \
      | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 -s \
      | samtools view | wc -l
83


It seems like I was wrong and COSMIC does store information about cell lines.

You need to download the files from Genotypes:

  Genotypes

  Files listing the SNP calls for each cell line identified by PICNIC
  analysis of Affymetrix SNP6.0 array data. Both a simple genotype (AA,
  BB - homozygous or AB - heterozygous) and a complex interpretation of
  the genotype are given (for example, in a triploid region of the
  genome the genotype maybe AAB). 'genotypes/.'

You need to register (and get a license) to be able to download the data. Presumably the cell lines are described in the heading of the file. Then you can sort them out by cell line.

Then once you know the germ line identifier look for the paper describing it to know more about their phenotype. (But maybe there is a database with phenotype of known cell lines)


Have a function that returns a single-end read (or pair) and another that does any relevant filtering on that. Wrap both of them in another function that checks if you've hit the end of the file. That's both simple, will handle any size of file, and allows single and paired-end reads.


How to select a power for a scale-free topology network

I would be concerned about the input data that produced these numbers. This dataset has a really odd shape with a dip at the start:

data.df <- read.table("soft_threshold.txt", header=TRUE);
png("scale-free-194.png");
plot(data.df$SFT.R.sq);
dummy <- dev.off();

Here's what I expect this to look like (from the supplementary information in this paper):

It's easier to select the power/threshold when the plot looks like the second situation. With the data as presented, it doesn't look like any value would be suitable. but if I were forced to pick, then I'd choose 12 because it is the first number above the magic threshold (ignoring the weirdness at the start). However, I'll admit that I don't understand the concepts behind the power, and am just following the recommended instructions, so can't explain why it is what it is.


Fasterq-dump: --split-spot or -concatenate-reads?

Fasterq comes from the latest version of sratools. So if you check the manual , it says the equivalence is:

fastq-dump SRRXXXXXX --split-3 --skip-technical
fasterq-dump SRRXXXXXX

In older versions of sratoolkit, if you use fastq-dump without specifying --split-3 for paired-end reads, you get the format mentioned, spotID.1 for forward, spotID.2 for reverse:
fastq-dump SRR6462984.sra
more SRR6462984.fastq

You get:
@SRR6462984.1 1 length=301
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGGATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG
+SRR6462984.1 1 length=301
#/AAAEE/EEEA/EAAEEAEE/EE/EEA/EEEEE/AAE/EE/EAEAE//EEAEE/EEEEE<EEEEE/EE<EE/AEEEE/EAEEEEA/E<EEAEAEEEEEEAEEEE<EEAEEEEEAAAEEE/EEEEE/E//AEEEE6EA66EE<<<EAAA/AAAAAE/E6E/EA/AE/AEEEEEEAAEEEEA6EEEEEEEEE/EEE/AEEAAAEEEAE/E/EEEEEA//E//EEE/EEE/A<//AE/EE//6A///AA/E//E/E/EE//AE//AAEEE////E/A<AAE/EEAE//A//6/A</<6/</E/
@SRR6462984.2 2 length=302
NTCAGCAGGCGCCGTAACTTCAGCAAAACTGTAGCTATAACTTTTCCGAACTTTATACCTAACTTCGCTTTCCTAGCACTGTATAAAATGCCTAAAAGAAATCAGCAGGCGCCGTCACTTCAAAAAATCAGTATGGATCAATCAGACGCGGAGCTTATTAAAACAGAGCTTATTAAAAATTATGTCTTTCTTCTTCTTTCAGCAATTCGCAACGTAAAATTCAAAATTTTGGACGCCGTAATTGCGCAGGAACCAATTCAAAGAGGCGGCGGCCAATTAATTCCAACGTCTCGAAGAGTGGC
+SRR6462984.2 2 length=302

With the latest version, it by defaults splits the paired end sra file into forward and reverse fastq (i.e --split-3):
fasterq-dump SRR6462984.sra
head -2 SRR6462984_*.fastq

And you get
==> SRR6462984_1.fastq <==
@SRR6462984.1 1 length=150
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGG

==> SRR6462984_2.fastq <==
@SRR6462984.1 1 length=151
ATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG

So if you need to get back the .1 and .2 format, use fastq-dump from the latest installation.
I would suggest sticking to the newest format, with 2 files, because fastq-dump will be deprecated in the future.


Difference in Calculation through GEO2R and difference in Article

Here are two example of calculations of logFC validated in GEO2R analysis which depends if values are already log2 transformed. During GEOQuery analysis there is calculation to transform where expression set is checked for value distribution and log2 transformation is applied-
ex <- exprs(gset)
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
          (qx[6]-qx[1] > 50 && qx[2] > 0)
if (LogC) { ex[which(ex <= 0)] <- NaN
  exprs(gset) <- log2(ex) }

Example 1:
=> Sample Values are already log2 Transformed
n(TumorSamples) = 45
n(ControlSamples) = 46

a = average(TumorSamples) = 9.199053111
b = average(ControlSamples) = 4.988139348

logFC = |a - b| = 4.210914

Example 2:
=> Sample Values are not log2 transformed
n(TumorSamples) = 2
n(ControlSamples) = 2

a = average(log2(TumorSample1) + log2(TumorSample2) ... ) = (log2(84) + log2(21.8)) / 2 = 5.419287

b = average(log2(ControlSample1) + log2(ControlSample2) ... ) = (log2(624) + log2(665)) / 2 = 9.331306

logFC = |a - b| = 3.91202


Filtering out all seqs with mutations of list2 from list1

Samtools/VCF is overkill for comparing amino acid sequences from a single protein. I would suggest making pairwise alignments between each sequence in list2 and the reference and making a list of mutations (2c and 6c in your example). Then repeat the procedure for list1, discarding sequences if they contain mutations in your stored list. This can all be done using biopython. Here's how to make the alignments to get you started:

from Bio import pairwise2
from Bio import SeqIO
seq1 = SeqIO.read("alpha.faa", "fasta")
seq2 = SeqIO.read("beta.faa", "fasta")
alignments = pairwise2.align.globalxx(seq1.seq, seq2.seq)


NCBI - edirect suite to download all genome sequences associated with a query - troubleshooting

The issue here was I was using an older version of eDirect. There is no command to check version or updates. Your best bet is if you experience problems uninstall and reinstall.


Downloading all COI sequences from BOLD fails

A straightforward way of getting the sequences from GenBank is from their website using this URL:
https://www.ncbi.nlm.nih.gov/nuccore?term=COI[Gene Name]

There you can click on "Send to" -> select "Gene Features" and "FASTA Nucleotide" -> Click "Create File".

There are >3 mio hits, so the download could take a while.


predict the foldability of single-stranded DNA molecules

You could also take a look at Infernal.

I think your best bet is to use RNA-focused programs; I'm not aware of a DNA-specific one, which makes some sense because DNA is typically double-stranded, so the secondary structure community is primarily focused on RNA. I would predict that secondary structure potential of single-stranded DNA would be similar to that of RNA.


What is the meaning of the Hmmer profile?

The tutorial in the documentation does indeed state what you quoted but a full description of the file format is given in Section 8 of the documentation, starting on page 106.
To summarise, as you noted, the number columns A-K is the alphabet size (the number of amino acids). The number of rows underneath is the number of nodes, which corresponds to the number of match states. The three numbers (-214, -210 and -18 in your example) are: the match emission, insert emission, and state transition values.
Caveat: The HMMER profile HMM format changes frequently between HMMER versions (even between minor versions) so there is no guarantee this information is correct for all the older file formats.

Quoting the relevant part of the documentation on page 109:

main model section
All the remaining fields are mandatory. The first
two lines in the main model section are atypical.‡ They contain
information for the core model’s BEGIN node. This is stored as model
node 0, and match state 0 is treated as the BEGIN state. The begin
state is mute, so there are no match emission probabilities. The first
line is the insert 0 emissions. The second line contains the
transitions from the begin state and insert state 0. These seven
numbers are: B → M1, B → I0, B → D1; I0 → M1, I0 → I0; then a 0.0 and
a ’’, because by convention, nonexistent transitions from the
nonexistent delete state 0 are set to log 1 = 0 and log 0 = −∞ = ‘’.
The remainder of the model has three lines per node, for M nodes
(where M is the number of match states, as given by the LENG line).
These three lines are (K is the alphabet size in residues):
Match emission line
The first field is the node number (1 . . . M). The
parser verifies this number as a consistency check (it expects the
nodes to come in order). The next K numbers for match emissions, one
per symbol, in alphabetic order. The next field is the MAP annotation
for this node. If MAP was yes in the header, then this is an integer,
representing the alignment column index for this match state
(1..alen); otherwise, this field is ‘-’. The next field is the CONS
consensus residue for this node. If CONS was yes in the header, then
this is a single character, representing the consensus residue
annotation for this match state; otherwise, this field is ‘-’. The
next field is the RF annotation for this node. If RF was yes in the
header, then this is a single character, representing the reference
annotation for this match state; otherwise, this field is ‘-’. The
next field is the MM mask value for this node. If MM was yes in the
header, then this is a single ’m’ character, indicating that the
position was identified as a masked position during model
construction; otherwise, this field is ‘-’. The next field is the CS
annotation for this node. If CS was yes, then this is a single
character, representing the consensus structure at this match state;
otherwise this field is ‘-’.
Insert emission line
The K fields on this
line are the insert emission scores, one per symbol, in alphabetic
order.
State transition line
The seven fields on this line are the
transitions for node k, in the order shown by the transition header
line: Mk → Mk+1, Ik, Dk+1; Ik → Mk+1, Ik; Dk → Mk+1, Dk+1. For
transitions from the final node M, match state M + 1 is interpreted as
the END state E, and there is no delete state M + 1; therefore the
final Mk → Dk+1 and Dk → Dk+1 transitions are always * (zero
probability), and the final Dk → Mk+1 transition is always 0.0
(probability 1.0).


TCGA: gene IDs not appearing and other concerns

I had to play around a little to understand the object, but here's what I got:
R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(TCGAbiolinks)
> query <- GDCquery(project = "TCGA-BRCA",
+                   data.category = "Gene expression",
+                   file.type  = "normalized_results",
+                   data.type = "Gene expression quantification",
+                   experimental.strategy = "RNA-Seq",
+                   barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
+                   legacy = TRUE)
--------------------------------------
o GDCquery: Searching in GDC database
--------------------------------------
Genome of reference: hg19
--------------------------------------------
oo Accessing GDC. This might take a while...
--------------------------------------------
ooo Project: TCGA-BRCA
--------------------
oo Filtering results
--------------------
ooo By experimental.strategy
ooo By data.type
ooo By file.type
ooo By barcode
----------------
oo Checking data
----------------
ooo Check if there are duplicated cases
ooo Check if there results for the query
-------------------
o Preparing output
-------------------

> GDCdownload(query, method = "api", files.per.chunk = 10)
Downloading data for project TCGA-BRCA
GDCdownload will download 3 files. A total of 1.306357 MB
Downloading chunk 1 of 1 (3 files, size = 1.306357 MB) as Wed_Aug_19_12_20_11_2020_0.tar.gz

Downloading: 580 kB

> data <- GDCprepare(query, summarizedExperiment = TRUE)
------------------
oo Reading 3 files
------------------
|====================================================|100%                      Completed after 0 s
------------------
oo Merging 3 files
------------------
Accessing grch37.ensembl.org to get gene information
Downloading genome information (try:0) Using: Human genes (GRCh37.p13)
Starting to add information to samples
 => Add clinical information to samples
 => Adding TCGA molecular information from marker papers
 => Information will have prefix 'paper_'
brca subtype information from:doi.org/10.1016/j.ccell.2018.03.014

> head(data)
class: RangedSummarizedExperiment
dim: 6 3
metadata(1): data_release
assays(1): normalized_count
rownames(6): A1BG A2M ... RP11-986E7.7 AADAC
rowData names(3): gene_id entrezgene ensembl_gene_id
colnames(3): TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07
  TCGA-AC-A23H-01A-11R-A157-07
colData names(83): barcode patient ... paper_PARADIGM Clusters
  paper_Pan-Gyn Clusters

I see from this representation that I can use metadata(data), assays(data), etc to get more information.
> head(SummarizedExperiment::assays(data))
List of length 1
names(1): normalized_count

Alright, so this is a named list. I know how to deal with lists: I can use either the index ([[1]]) or the name ($normalized_count)
> head(SummarizedExperiment::assays(data)$normalized_count)
             TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-01A-11R-A157-07
A1BG                             121.9879                      38.5795                      13.6874
A2M                             9703.8622                   37640.8373                    4767.0490
NAT1                            8488.0727                     114.8582                     270.4753
NAT2                               6.0583                       4.2938                       2.3804
RP11-986E7.7                   33016.6641                   11308.3460                    3894.9639
AADAC                              1.1359                     219.6976                       1.1902

I can now get more information on the samples using colData(). I'm truncating the output below because there are 83 columns.
> SummarizedExperiment::colData(data)
DataFrame with 3 rows and 83 columns
                                                  barcode      patient           sample shortLetterCode          definition sample_submitter_id         sample_type       state
                                              <character>  <character>      <character>     <character>         <character>         <character>         <character> <character>
TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046 TCGA-AN-A046-01A              TP Primary solid Tumor    TCGA-AN-A046-01A       Primary Tumor    released
TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-11A              NT Solid Tissue Normal    TCGA-AC-A23H-11A Solid Tissue Normal    released
TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-01A              TP Primary solid Tumor    TCGA-AC-A23H-01A       Primary Tumor    released
                                                        sample_id   is_ffpe sample_type_id  tissue_type submitter_id days_to_collection oct_embedded                pathology_report_uuid
                                                      <character> <logical>    <character>  <character>  <character>          <integer>  <character>                          <character>
TCGA-AN-A046-01A-21R-A034-07 57aff09f-0e97-4e65-a227-3dc7a8516367     FALSE             01 Not Reported TCGA-AN-A046                 27         true 1304FB17-A20A-4EBC-9CC8-1554808AC1F6
TCGA-AC-A23H-11A-12R-A157-07 7df59ca8-7e51-4581-9d7f-8bba0395ce17     FALSE             11 Not Reported TCGA-AC-A23H                478        false                                   NA
TCGA-AC-A23H-01A-11R-A157-07 d7e3b628-d5fd-4e79-9c4a-6409330fb8a7     FALSE             01 Not Reported TCGA-AC-A23H                478        false A7C7D409-D086-4A9B-8C8F-E7E231D5891D

I can get more information on the genes using rowData():
> SummarizedExperiment::rowData(data)
DataFrame with 19947 rows and 3 columns
                  gene_id entrezgene ensembl_gene_id
              <character>  <integer>     <character>
A1BG                 A1BG          1 ENSG00000121410
A2M                   A2M          2 ENSG00000175899
NAT1                 NAT1          9 ENSG00000171428
NAT2                 NAT2         10 ENSG00000156006
RP11-986E7.7 RP11-986E7.7         12 ENSG00000273259
...                   ...        ...             ...
RASAL2-AS1     RASAL2-AS1  100302401 ENSG00000224687
LINC00882       LINC00882  100302640 ENSG00000242759
FTX                   FTX  100302692 ENSG00000230590
TICAM2             TICAM2  100302736 ENSG00000243414
SLC25A5-AS1   SLC25A5-AS1  100303728 ENSG00000224281

You want the assay data as a data.matrix, I guess. That'd be:
> data_as_matrix <- data.matrix(SummarizedExperiment::assays(data)$normalized_count)

My session Info:
> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] TCGAbiolinks_2.16.1

loaded via a namespace (and not attached):
  [1] bitops_1.0-6                matrixStats_0.56.0          bit64_0.9-7.1               doParallel_1.0.15           RColorBrewer_1.1-2          progress_1.2.2
  [7] httr_1.4.2                  GenomeInfoDb_1.24.2         backports_1.1.8             tools_4.0.2                 R6_2.4.1                    DBI_1.1.0
 [13] BiocGenerics_0.34.0         colorspace_1.4-1            tidyselect_1.1.0            gridExtra_2.3               prettyunits_1.1.1           bit_1.1-15.2
 [19] curl_4.3                    compiler_4.0.2              rvest_0.3.6                 Biobase_2.48.0              xml2_1.3.2                  DelayedArray_0.14.0
 [25] rtracklayer_1.48.0          scales_1.1.1                survMisc_0.5.5              readr_1.3.1                 genefilter_1.70.0           askpass_1.1
 [31] rappdirs_0.3.1              stringr_1.4.0               digest_0.6.25               Rsamtools_2.4.0             foreign_0.8-80              R.utils_2.9.2
 [37] rio_0.5.16                  XVector_0.28.0              pkgconfig_2.0.3             dbplyr_1.4.4                readxl_1.3.1                rlang_0.4.7
 [43] ggthemes_4.2.0              RSQLite_2.2.0               generics_0.0.2              zoo_1.8-8                   jsonlite_1.7.0              BiocParallel_1.22.0
 [49] zip_2.0.4                   car_3.0-8                   dplyr_1.0.0                 R.oo_1.23.0                 RCurl_1.98-1.2              magrittr_1.5
 [55] GenomeInfoDbData_1.2.3      Matrix_1.2-18               Rcpp_1.0.5                  munsell_0.5.0               S4Vectors_0.26.1            abind_1.4-5
 [61] lifecycle_0.2.0             R.methodsS3_1.8.0           stringi_1.4.6               carData_3.0-4               SummarizedExperiment_1.18.1 zlibbioc_1.34.0
 [67] plyr_1.8.6                  BiocFileCache_1.12.0        grid_4.0.2                  blob_1.2.1                  parallel_4.0.2              ggrepel_0.8.2
 [73] forcats_0.5.0               crayon_1.3.4                survminer_0.4.8             lattice_0.20-41             haven_2.3.1                 Biostrings_2.56.0
 [79] splines_4.0.2               GenomicFeatures_1.40.0      annotate_1.66.0             hms_0.5.3                   knitr_1.29                  pillar_1.4.6
 [85] ggpubr_0.4.0                GenomicRanges_1.40.0        ggsignif_0.6.0              codetools_0.2-16            biomaRt_2.44.1              stats4_4.0.2
 [91] XML_3.99-0.5                glue_1.4.1                  downloader_0.4              data.table_1.13.0           vctrs_0.3.2                 selectr_0.4-2
 [97] foreach_1.5.0               cellranger_1.1.0            gtable_0.3.0                openssl_1.4.2               purrr_0.3.4                 tidyr_1.1.0
[103] km.ci_0.5-2                 assertthat_0.2.1            ggplot2_3.3.2               openxlsx_4.1.5              xfun_0.16                   xtable_1.8-4
[109] broom_0.7.0                 rstatix_0.6.0               survival_3.1-12             tibble_3.0.3                iterators_1.0.12            KMsurv_0.1-5
[115] GenomicAlignments_1.24.0    AnnotationDbi_1.50.1        memoise_1.1.0               IRanges_2.22.2              ellipsis_0.3.1


Why does IGV turn a VCF QUAL of "." into -10?

Looks like they are treating the "." as a value of 1.  According to the VCF spec the quality is "−10log10 prob".  Where are you getting your QUALs from?  I have never seen the quals as a "." and IGV probably needs to do some translation into an int but a QUAL of -10 makes no sense.

I would report it and let them explain.


Converting FASTA Sequences to SNP FASTA Format

Your question is not very clear. Where does the flanking sequence in :
>rs11540652
GAACC[A/G]GAGGCCC

come from? I'll assume its meant to be the reference sequence and that speciesA is the reference and that the first nuc in the [] is from the reference... I don't know an existing existing program but here's a bespoke one.
speciesA='TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
speciesB='TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
counter=0
for i, nuc in enumerate(speciesA):
    if nuc != speciesB[i]:
        counter+=1
        print ('>rs'+str(counter))
        print (f'{speciesA[i-5:i]}[{nuc}/{speciesB[i]}]{speciesA[i+1:i+6]}')

Gives:
>rs1
AAAGT[G/A]ATGTA
>rs2
TGTAA[T/C]CCTAC
>rs3
CAAAA[G/-]TTATA
>rs4
TTATA[A/G]GAACT
>rs5
CATTA[T/C]ATCAT
>rs6
TTATA[T/C]CATCA
>rs7
AACCA[C/T]TAAAC
>rs8
CCACT[A/G]AACTT
>rs9
CTAAA[C/T]TTCAA
>rs10
AACTT[C/T]AAACT
>rs11
TCAAA[C/T]TCCAT
...

As for 'I presume that I need assembled mitochondrial reference genomes for each species - anything else?' These sequences have clearly been aligned. So there's that.
I hope that helps.


How to concatenate FASTAs with respect to chromosome?

Using the TblToFasta and FastaToTbl scripts I have posted before, combined with classic *nix utils, you can do:

$ join -a1 -a2 <(FastaToTbl dna1.fa | sort) <(FastaToTbl dna2.fa | sort) |
    sed 's/ /=/2' | TblToFasta | tr = '\n'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation

The FastaToTbl will put everything on the same line, with the sequence's ID at the beginning of the line, then a tab and then the sequence (no matter how many lines of sequence you have). So, we do this to both files, sort them and pass that as input to join, which will join its input on lines where the 1st field is the same (by default). The -a1 -a2 means "also print non-matching lines from both files" so we don't miss cases where a sequence is in only one of the two files.

Since the join will add a space between the two joined lines, we deal with that by replacing the second space on each line with a = (sed 's/ /=/2') and then pass the whole thing to TblToFasta to revert to fasta format. Finally, we use tr to replace the = with a newline so you can see where each sequence came from as you show in your example.

For longer sequences where you want to avoid sorting the input, you can do:

$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2
         }
         END{
              for(chr in seq){
                printf ">%s\n%s\n", chr, seq[chr]
              }
         }'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation

Here, the tbl format sequences (see above) are passed to an awk script which does the work:

seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2 : if the 1st field ($1, the sequence identifier) is already in the array seq, then set the value of seq[$1] to the current value (so whatever sequence has been already seen for this ID), a newline and the 2nd field ($2, the sequence). If it hasn't been seen before, then set the value of seq[$1] to the 2nd field, the sequence.
END{ for(chr in seq){ printf ">%s\n%s\n", chr, seq[chr]} : after processing everything, go through each element in the array seq and print out the data in fasta format.

However, that will mean keeping both files in memory which may be a problem for huge sequences. You can avoid that by using temp files:

$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            if(!seen[$1]++){
                print ">"$1 > $1".tmp"
            }
            print $2 >> $1".tmp"
        }' && cat *.tmp > out.fa && rm *.tmp
$ cat out.fa
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation

This is the same basic idea as the script above, only instead of storing the information in an array, we print it as it comes into temp files.

if(!seen[$1]++){print ">"$1 > $1".tmp"} : if this is the first time we've seen this ID, print it as a fasta header into the file $1".tmp" (e.g. chr1.tmp).
print $2 >> $1".tmp" append the current sequence to the corresponding temp file.

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (rm *tmp).

As for why you would want to do this, an obvious example is reconstructing a reference genome from short sequences. You would want all the sequences of the same schromosome together. Of course, you would also need to assemble them in the right order, which you aren't doing here.

Alternatively, you might want to concatenate them this way so as to avoid having the extra, unneeded lines (>chrN) which would make your file larger. Or you would want to be able to retrieve all sequences with the same ID at once.

But I don't really know why you would want to do it. You're the one asking for it, after all! :)


How to perform DE analysis for each sample

Treat your disease samples as individual groups and then follow the normal routine in limma to use contrasts to compare two group (the control group versus individual disease samples). Note that you need to think long and hard about what these results then mean. Normally we compare groups because then the results should generalize to other samples. That will not be the case when you look at individual disease samples, so I suggest you do some clustering of the samples to see if you perhaps have any coherent clustering of the disease samples into subsets that could then be used as more useful groups in DE testing.


Methylation can be cell-specific, which makes it difficult to evaluate accuracy on a bulk-cell level (even within the same tissue). How can you tell that the differences you're seeing are due to platform differences, or due to biological variation?
I find that adding more haystacks doesn't help much in working out the truth of a dataset. If you want to investigate biological truth, it would be better to create a biological system with a known methylation (or demethylation) pattern and test that.