How to access organism IDS, taxon nodes from UniProt?

These are NCBI taxonomy ID's rather than anything specific to UniProt. One way to parse them would be using BioPython via Bio.Entrez, for example:

from Bio import Entrez

def get_taxon_name_from_NCBI_id(id_):
    Entrez.email = 'name@provider.com'
    handle = Entrez.efetch(db='taxonomy', id= id_, retmode='xml')
    records = Entrez.read(handle)[0]
    return records['ScientificName'], records['TaxId']

print(get_taxon_name_from_NCBI_id(9606))
# ('Homo sapiens', '9606')

Of course you could query Entrez directly:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy&id=9606

Or even more basically do a manual search of NCBI taxonomy:
https://www.ncbi.nlm.nih.gov/taxonomy/?term=9606


Albacore basecalling running but outputs 0 reads

I believe you have multiread fast5 files, as generated by recent versions of MinKNOW. A conversion script is available to convert these to the older, one-read one-file format. AFAIK Albacore is being deprecated and no version will be made supporting the multiread fast5 files, as Albacore is being replaced by Guppy.

The command line to convert from multi read fast5s to single read fast5s is as follows:

$multi_to_single_fast5 -i /path/to/fast5/ -s /path/to/save/output/


How to compute LD for pairs of variants with Plink

Yes, you can do this with plink.

You can use the flag --ld-snp <variant_id> to calculate LD between your SNP of interest and all other SNPs.
If you have multiple SNPs of interest, and want each of those compared with every SNP in your dataset, use --ld-snps <snp1-snp100,snp102,snp104>
A final option to calculate LD between a specific pair of SNPs is --ld <variant_id1> <variant_id2>

This the relevant plink manual page


Important genes beyond PAM50 for breast cancer classification

There is a CSV table in this paper with 33 sets of genes.


Alignment for predicting DNA hybridization?

As a rough and ready and very easy to perform test I would put the entire universe of dsDNA through Primer3 (details below). This could be done by "splicing" your probe of interest onto each "chunk" of dsDNA you pass through primer3 and restrict the primer search to the 20bp primer you are interested in. In other words you are forcing the program to consider the internal and trans-hybridisation of every possible dsDNA against your probe. If you probe is not recovered from a given search you would flag this on a subsequent parsing script and then that "chunk" of DNA would be assessed further. All the components required are in the link below. You would need to provide your own wrapper, but its fairly easy.

Primer3 will assess the internal thermodynamics of your probe and it has loads of options to assess the full range of behaviour, for example 3' hybridisation (which could be extended because you want 5' to 3'). It would be easy to construct a parameter grid around all the parameters of interest. Ultimately it would provide an additional layer of certainty, and everyone knows and trusts this algorithm (its also written in C, so its fast).

The description of all the components is here , my specific purpose for customising the pipeline is different.

Bowtie2 will take a lot of RAM, I am not sure I would use it for this purpose, particularly if RAM is limiting (its used for genome assemblies against a reference). I'd ask the users whether Blat would be useful, because that is friendlier to the calculation you are performing. Anyway there some heavy duty Bowtie2 users around and their input would be useful.


This was addressed by the Seurat developers here:

  if you have TPM counts, I suggest you don't use
  Seurat::NormalizeData(), since TPM counts are already normalized for
  sequencing depth and transcript/gene length. Note that
  Seurat::NormalizeData() normalizes the data for sequencing depth, and
  then transforms it to log space. If you have TPM data, you can simply
  manually log transform the gene expression matrix in the object@data
  slot before scaling the data.

And here:

  If you want to run your own normalization function, please place the
  data in log scale before placing it into object@data

And here:

  If you're supplying already normalized values to CreateSeuratObject,
  you should skip the normalization step. You will likely see a warning
  when you run ScaleData (which is checking whether you normalized using
  Seurat). You can ignore this or turn it off by setting check.for.norm
  to FALSE.

This is a lot of quotes, but I found them all to be helpful in different ways.


SNPEFF genome build with specific codon usage

codon table issue

I solved the issue with codon tables, the correct way to specify a table is the following:

Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : Bacterial_and_Plant_Plastid

Note that one should remove the leading prefix 'codon.' from the name of the tables chosen from the config file

Concerning accessory files:

data to reproduce this issue:

baseurl="ftp://ftp.ensemblgenomes.org/pub/plants/release-46"
wget ${baseurl}/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ${baseurl}/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gff3.gz
wget ${baseurl}/variation/vcf/arabidopsis_thaliana/arabidopsis_thaliana_incl_consequences.vcf.gz

I suspect that the format used by Ensembl is not 100% OK for SnpEff; I therefore simplified the fasta headers to keep only the first argument (name) and remove the remaining of the line

eg.

>AT3G05780.1 pep chromosome:TAIR10:3:1714941:1719608:-1 gene:AT3G05780 transcript:AT3G05780.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:LON3 description:Lon protease homolog 3, mitochondrial [Source:UniProtKB/Swiss-Prot;Acc:Q9M9L8]

becomes

>AT3G05780.1

Despite this modification, I still have an issue with

java -jar /opt/biotools/snpEff/SnpEff.jar build -gff3 -v Arabidopsis_thaliana_TAIR10_46 2>&1 | tee -a log.txt

a LOT of warning about not being able to add UTR while it seems that UTR's are well present in the GFF3 file

  Reading GFF3 data file  :
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  WARNING: Cannot find transcript 'transcript:AT1G01010.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  18    '1  araport11   five_prime_UTR  3631    3759    .   +   .   Parent=transcript:AT1G01010.1'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  18    '1  araport11   five_prime_UTR  3631    3759    .   +   .   Parent=transcript:AT1G01010.1'
  WARNING: Cannot find transcript 'transcript:AT1G01020.3'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  69    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.3'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  69    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.3'
  WARNING: Cannot find transcript 'transcript:AT1G01020.5'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  88    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.5'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  88    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.5'
  WARNING: Cannot find transcript 'transcript:AT1G01020.4'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  109   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.4'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  109   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.4'
  WARNING: Cannot find transcript 'transcript:AT1G01020.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  128   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.1'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  128   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.1'
  WARNING: Cannot find transcript 'transcript:AT1G01030.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line 155  '1  araport11   three_prime_UTR 11649   11863   .   -   .   Parent

nothing good coming out of the cds.fa input although I treated it with the same name fix as the pep file

  00:00:19  CDS check (FASTA file):
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/cds.fa'

  00:00:19  Checking database using CDS sequences 00:00:19  Reading CDSs
  from file
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/cds.fa'...
  00:00:20  done (75976 CDSs). 00:00:20 Comparing CDS...    Labels:         '+' :
  OK        '.' : Missing       '*' : Error
    ....................................................................................................
    ....................................................................................................
    ....................................................................................................

  many more lines with dots
    ....................................................................................................
    ....................................................................................................
    ....................................................................................................
    ...............................     CDS
  check:    Arabidopsis_thaliana_TAIR10_46  OK: 0   Warnings: 0 Not found:
  104147    Errors: 0   Error percentage: NaN% Fatal error: No CDS checked.
  This is might be caused by differences in FASTA file transcript IDs
  respect to database's transcript's IDs. Transcript IDs from database
  (sample):     'transcript:AT3G12250.7'    'transcript:AT5G18570.1'
    'AT3G27810.1'   'AT3G27810.2'   'transcript:AT2G16810.1'
    'transcript:AT3G01170.1'    'AT4G37030.1'   'transcript:AT1G79940.4'
    'transcript:AT2G44850.1'    'AT1G49540.1'   'AT1G49540.2'
    'transcript:AT4G35290.5'    'transcript:AT4G15840.2'
    'transcript:AT4G27500.1'    'AT4G37022.1'   'transcript:AT2G40008.1'
    'AT3G27820.1'   'transcript:AT4G24265.1'    'transcript:AT2G20618.1'
    'transcript:AT4G18375.2'    'AT4G37020.2'   'AT4G37020.1' Transcript IDs
  from database (fasta file):   'AT4G02541'     'AT1G12410.1'   'AT4G26510'
    'AT4G02540'     'AT3G04500.1'   'AT3G04500.2'   'AT3G04500.3'   'AT2G36250'
    'AT5G63350'     'AT2G36255'     'AT2G12280'     'AT1G03970.1'   'AT5G18390.1'
    'AT3G04490.4'   'AT5G38386'     'AT3G04490.2'   'AT3G04490.3'
    'AT3G04490.1'   'AT5G38380'     'AT4G02530'     'AT4G26520'     'AT1G69280.2'

link to the full stdout log after running the command without cds.fa

https://www.dropbox.com/s/gxa4e533p64wxm9/SnpEff_build-log.txt?dl=1


Filtering the expression values of a gene from GEO

You're looking for the Columns() function:

> Columns(gds3782)
     sample        disease.state           description
1  GSM524151 non-diabetic control  Value for GSM524151: beta-cells_non-diabetic condition_donor1; src: human beta-cells_non-diabetic
2  GSM524152 non-diabetic control  Value for GSM524152: beta-cells_non-diabetic condition_donor2; src: human beta-cells_non-diabetic

Note that it's probably easier to work with a more common structure like and expressionSet:

eset = GDS2eSet(gds3782)

You'd then use the pData() function to return the same thing.


How to achieve blast results according to the intuitive interpretation of `-max_target_seqs`?

The -max_target_seqs is also related to -evalue which does not work as one would think. Looking at the blast news one finds that since version BLAST 2.2.27+ released, one should use the -max_target_seqs:

  4.) The output for reports without separate descriptions and alignments sections (all –outfmt greater than 4) should now use –max_target_seqs to control the output rather than –num_descriptions and –num_alignments.

So it seems that num_alignments is also misunderstood. I came to the conclusion that one cannot expect to get the best n alignments, unless you get all and filter yourself the best ones.


As always with bioinformatics, any reasonable answer will be dependendant on what you want to do with this sequence, and you haven't provided enough information in that regard.

For most exploratory analysis BLAST is good enough. You don't have to follow any existing methods. The benefit of using NCBI's web BLAST server is that it provides you with a very quick way to find a match from the sequence information in GenBank.

BLAST does a local alignment. If you have the full 16s sequence, and it's a novel organism, you might get different match results from a program that does global alignment, but then again you might get different results from a different local alignment tool.


Differences between NetAffx Hg-U133 Plus 2.0 Annotation file versions

After checking the files myself, it does seem like gene identifiers indeed have been updated. Some that were missing have been added with identifiers and some have additional identifiers now.

Among some other things that were updated from 2007 include:

Genome version NCBI Version 36.1 => GRCh37
Alignments
Gene title
Gene symbol

Thought this would be a good reference for those who may be looking into older researches involving microarrays. I will be sure to check the files more carefully before posting questions like this one.


A basic question about phylogenetic inference with Bayes-MCMC

In short, the posterior is a probability value between 0 and 1 which is calculated from Bayes rule, which In phylogenetic usage is given as the product of the prior probability and the likelihood divided by a normalising constant. The posterior probability allows you to state the probability of some hypothesis given the data at hand (as opposed to the likelihood parlance, which is probability of data given hypothesis).

In a phylogenetic setting, the sampler works (roughly) like this: propose some value of a parameter (topology, branch lengths, etc.) and use the specified prior and calculated likelihood for that parameter value given the data to calculate the unnormalised posterior. Next, propose a new value, and repeat the same step. Compare the two values, and if the posterior probability of the second is greater, accept this parameter value and repeat the cycle again. In this way the Markov Chain moves through posterior space until it reaches a reasonably stable optimal approximation of the posterior distribution. Here, we want to get as many samples as possible to maximise the accuracy of our posterior. The maximum a posteriori (MAP) value of the parameter and its credible intervals, etc, will come from the distribution of these samples.


gmx mdrun error: "Error in user input: Invalid input values In option s"

If the command line says "command not found", it means that the binary is not executable or not found in any of the paths mentioned in the $PATH environmental variable. You can check for example this tutorial for details on how it works and how to resolve problems related to this variable.
Which probably means, that you need to go a step back, to the installation instructions of gromacs and try to install the mdrun binary. However, if the mdrun was not installed successfully, it's very likely that more stuff of the package will be missing, so it might be worth reinstalling it all.
P.S. It's not obvious to me if getting mdrun it will resolve the first posted problem (I have never worked with gromacs). If you would share the link to the research you have done, that would help others to help you.


How is consensus alignment for OLC assembly usually implemented?

For short reads, mainstream assemblers don't have a consensus step. The contig path directly spells the consensus sequence.

For long reads, the most popular way to derive a consensus is via partial-order alignment (POA). This is used by error correctors (similar to consensus) like pbdagcon and falcon_sense, and consensus tools like canu, racon, wtdbg2 and nanopolish. Basically, you start with an empty directed acyclic graph (DAG), align each read in turn to the graph and augment the DAG by inserting unaligned segments. When you align all reads onto the DAG, the most weighted path is taken as the consensus. You should read Jared Simpson's blog on this topic. It is much more detailed than "basically". You should also read the papers on all the mentioned tools above if you want to implement by yourself.

An alternative approach is Gene Myers' realigner (PDF here). It is not often used, though.


How to retrieve SNPs data of different humans?

If you mean that you want snps from individuals, instead of all together, you can find them in 1000 genomes. Here different individuals from different populations are sequenced and variants are called, including snps.

Maybe you should add to your question what your goal is, what would you like to find?


What is the correct way to map Hi-C data with bwa mem?

See here, in particular slide #10 from the tutorial:

bwa mem -5SPM ref.fa read1.fq read2.fq > out.sam

Here, -SP disables pairing. The Aiden lab and this paper also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.


Whenever you are wondering about things like this, just look up the identifier on the Ensembl web page. If you look up ENSMUSG00000083840, you will see:

  This identifier is not in the current EnsEMBL database

In this case, it doesn't seem to have been replaced by anything and it wasn't a protein coding gene, so it looks like it was simply a wrong annotation which has been removed.


Is there a computational method/method/way to predict if a cell population in a tissue immigrated or is enriched locally?

This may not be a viable answer for all experimental contexts, but you can compare expression profiles of individual cells against datasets that contain known reference cell types (migratory and not). We have recently published a web-based tool called CIPR (cluster identity predictor) that helps with cluster annotations in scRNAseq experiments. There are various reference datasets available, but you can also provide a custom reference by preparing data from databases like GEO.
You can read the manuscript for more information about how the algorithm works, but in a nutshell, it compares differentially expressed genes or global expression patterns within your data to known reference datasets and creates quick visual outputs. If your reference contains both migratory and non-migratory cell types, maybe CIPR can help differentiate them. For other identity prediction tools, check out this article is a good resource that compares different cluster calling algorithms.


How to align output of grep --color=always? (To QC fasta/fastq files)

Perhaps, grep is not the best tool to use in this case, but it should be in principle possible by using grep & sed. Here is an example showing three symbols around a match.

zcat My_Hiseq_Data.fq.gz | \
      grep -Eo '.{0,3}GATCGATC.*' | \
      sed -En 's/.*/   \0/; s/.*(.{3}GATCGATC.{0,3}).*/\1/p' | \
      grep --color=always GATCGATC

Here are some explanations:

zcat will just decompress a file.
grep will output all the matches with 3-symbol context (or less). This is needed because a line can have multiple matches.
first sed command will add 3 symbols before all the matches.
second sed command will output the sequence with 3 surrounding characters on the left and max 3 characters on the right. This is an equivalent of grep -Eo '.{3}GATCGATC.{0,3}'.
Finally, now that all matches are aligned, grep will color your sequence.

This code is possible to adapt to include more symbols of context.

I would recommend doing this in another language, though. E.g., in Python:

#!/usr/bin/env python3
import sys
import re

HL = '\033[93m'
ENDC = '\033[0m'

filename = sys.argv[1]
seq = sys.argv[2]
context = 10

seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (context, sys.argv[2], context))

for line in open(filename):
    line = line.rstrip()
    for m in seqre.finditer(line):
        print(m.group(1).rjust(context),
              HL,
              m.group(2),
              ENDC,
              m.group(3),
              sep='')

UPD

A more advanced Python version inspired by this answer with support for multiple files, gzip, etc.

#!/usr/bin/env python3
import sys
import re
import argparse
import gzip
from termcolor import  colored

def align_matches(f, seqre, context):
    # for every line
    for line in f:
        # get rid of newline character
        line = line.rstrip()
        for m in seqre.finditer(line):
            # for every match print right-aligned left-context,
            # match colored in red, and right-context
            print(m.group(1).rjust(context),
                  colored(m.group(2), 'red'),
                  m.group(3),
                  sep='')

# this function will check file extension to be able
# to open text files as well as gzip-compressed files
def myopen(fn):
    if fn.endswith('.gz'):
        return gzip.open(fn, 'rt')
    else:
        return open(fn, 'rt')

if __name__ == '__main__':
    # define script arguments
    parser = argparse.ArgumentParser(description='Aligned grep')
    parser.add_argument('pattern', type=str,
                        help='a valid regular expression')
    parser.add_argument('input', type=myopen, nargs='+',
                        help='files names, could be a gzip compressed file ending with .gz')
    parser.add_argument('--context', '-c', default=30, type=int,
                        help='number of context characters to show')
    args = parser.parse_args()

    # compile a regular expression with a match surrounded
    # by at most args.context number of characters
    seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (
        args.context, args.pattern, args.context))

    # for every input file run align_matches
    for f in args.input:
        align_matches(f, seqre, args.context)


BWA: Detecting Variation between Reference Genome and Short-Read Sequences

I haven't done any whole-genome STR analysis from NGS data myself, but are aware of others that have used lobSTR for this. There's also a recent paper [here] that compares a few different STR analysis packages (i.e. RepeatSeq, LobSTR, HipSTR, GangSTR). Here's the concluding paragraph:

In conclusion, all these tools are built to genotype STRs but
have different strengths and weaknesses. Based on our analysis there is
no clear overall winner. RepeatSeq and HipSTR are the best when
considering genotyping error rate even with low coverage. On the other
hand, GangSTR has an advantage because it is the only tool among them
that can call alleles longer than the read length but shows higher error
rate unless looking at only the enclosed class of reads, which in turn
would lose the GangSTR's advantage of picking up long genotypes. In
addition, GangSTR is the newest tool and so comes with reference files for
different reference builds that are periodically updated according to
the tool's webpage. The correct choice of a tool and the
subsequent filtering depends on the aim of the analysis, and might be
influenced by available hardware resources and time limit for running
tools.

In other words, without more information about the specific problem you're trying to solve, it's not clear what the best tool to use is (this is a common issue with Bioinformatics problems).


You can extract fasta sequence from a blastdb constructed from a fasta file using blastdbcmd which should be installed when you install blast/makeblastdb.

blastdbcmd -entry all -db <database label> -out <outfile>

If you had a database called my_database which contained the files:

my_database.nhr
my_database.nsq
my_database.nin

and you wanted your fasta output file to be called reference.fasta you would run the following:

blastdbcmd -entry all -db my_database -out reference.fasta


Mark BWA-SW split alignments in output for long reads

From BWA manual, it is advised to use BWA-MEM in your longs read case (unless you have a good reason not to):

  BWA-MEM and BWA-SW share similar features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and more accurate.

I guess that BWA-SW returns several lines per transcript (based on local alignments). In that case, the issue has already been reported. The user asks:

  How can I differentiate between a read that has multiple hits in the reference and a read that has a split alignment?

The author replies:

  In principle, it is hard to tell. Bwa-mem makes an arbitrary choice.
  Mappings flagged with 0x100 are multiple hits. Mappings flagged with
  0x800 are from a chimeric alignment.

So you could actually filter by flag. But I am afraid the output could include something you do not necessarily want to (since there is no specific flag for the split reads).

Maybe a better solution: you could switch to some mapper that will report the proper split in the hit of the read map, such as minimap2 from the same author. It is suitable for long reads and comes with a splice -x feature.

Then, the laziest solution I see here would be to remove the split alignments exploiting the sam format and in particular the CIGAR string. Meaning, to remove the sam fields containing a N in their CIGAR (aka, skipped region from the reference). Something like this would work:

cat alignment.sam | awk '$6 !~ /N/ {print $0}'

Hope this can be somewhat helpful.


Pipeline (Pypiper) returns error, but individual command works fine

Slurm is killing the process because it's using more memory per core than the default allocation (that's what slurmstepd: error: Exceeded step memory limit at some point means). I've never used pyPiper (like seemingly most people, I use snakemake), but presuming that you have it running either srun or sbatch at some point then specify something like --mem-per-cpu 10G.


How to remove zero value of gene on FeatureScatter plot using Seurat?

It will not improve your correlation most likely, unless the zeros are not missing at random:

library(Seurat)
data <- pbmc_small
data[["percent.mt"]] <- PercentageFeatureSet(data, pattern = "^MT-")
data <- NormalizeData(data)
genes <-c("HLA-DRA","LYZ")
countdata = pbmc_small[["RNA"]]@data
countdata <- t(as.matrix(countdata)[genes,])
plot(countdata)

Now you need to subset the matrix, and remove the rows where at least one of the entries is 0, you can use rowSums(..==0)==0, meaning no zeros in the rows:

plot(countdata[rowSums(countdata==0)==0,])

cor(countdata[rowSums(countdata==0)==0,])
          HLA-DRA       LYZ
HLA-DRA 1.0000000 0.2338371
LYZ     0.2338371 1.0000000


No. A perfect phylogeny is one such that all characters evolve on the tree with no homoplasy, i.e. for a binary character, changes occur once from 0 -> 1 but never from 1 -> 0. A maximum parsimony phylogeny may produce a perfect phylogeny, but typically for real datasets some degree of homoplasy is required to explain character patterns.


Prevent Snakemake from changing wildcards: sample for every rule

Placing two wildcards directly behind each other will make it impossible for snakemake to determine where one wildcard ends and the other begins. So it is simply taking the longest possible match for the first wildcard (sample) and then one character for the second (read). This will also cause the _ to be included in the wildcard value sample.

I'm assuming that this same pattern is used for the inputs of another rule to refer to the trim_galore_pe output. Hence why the _ is included in the wildcards there.

You might be able to fix this, simply by adding a _ in-between {sample} and {readgroup} in your first rule (and any other rules using this same pattern):

rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}_{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"


What samples can be used for a Mutect2 Panel of Normals

As suggested by @burger in their comment, you can indeed use solid tissue samples as normals for your PON. You are clearly aware, but I want to spell this out for others, that there are a few caveats that go with the practice of using solid tissue normals rather than blood normals. Mainly, if your normals are "tissue-adjacent", ie sampled from an area of tissue adjacent to the tumor, there is a higher risk of including contaminating tumor cells compared to blood normals. As a result you might have real somatic variants that end up filtered because they are seen in the normal. Though that typically acts more at the level of the direct Tumor-Normal comparison than at the level of the PON, since the PON draws from many normal samples and the chance of seeing the same contaminating somatic mutation in many such samples should be quite low. So ultimately that should be fine.

It's a very common practice to use whatever matched normals you have available to build your PON, and there doesn't seem to be any drawback to doing that as long as you have enough samples (at least a few dozen, ideally a few hundred or more). If you don't have enough samples in your own study, then you should seek to use other samples from young and healthy people that have been sequenced using the same technology (same library prep kits/protocols, same chemistry version etc) since the goal of the PON is to weed out recurrent technical artifacts.


How to interpret methylation calls from Bismark on opposite strands?

The strand the bismark reports is related to the strand from which the read originated, not necessarily how it's aligned. So, alignments on the + strand shouldn't have calls overlapping those on the - strand, since you can't have a C in the same place on the same strand. One should often see Z/z next to each other on opposite strand, like in your example, since these are CpG (so it should be a Z/z on the + strand and then a Z/z on the following base on the - strand). Thus, in your example you have two reads supporting unmethylation for the CpG as a whole (1 read for each of the Cs).

The most confusing thing about BSseq is that reads (or paired-end reads) are only ever informative for a single strand. This is actually true for all sequencing with Illumina instruments (single-stranded fragments are loaded, after all), but since strands are almost always complementary we can usually ignore this fact.

As an aside, you might find MethylDackel useful (full disclosure, I wrote it). It'll be much faster at extracting methylation calls than bismark and supports nice things like excluding regions of bias methylation and likely variant positions.


From looking at this publication "Database indexing for production MegaBLAST searches" which although deals with MegaBLAST and it's indexing method, it would appear that this was rolled into the NCBI C++ toolkit as the command makembindex so I'm inclined to think the current implementation must have been derived from this strategy. Looking at the relevant part of the current NCBI C++ toolkit I find no reference to makembindex but only makeblastdb which I assume has superseded the older command: blastdb Directory Reference

Ulimtately the true details of the index generation lie in the C++ source code for makeblastdb can be found here.  However following on from that I think the further details of the index generation may be in one of the following files which are includes for the main .cpp file:

seqdb.hpp

build_db.hpp

seqdbcommon.hpp

writedb.hpp

It would appear from seqdb.hpp that:

The index files have extension .pin or .nin and sequence files .psq or .nsq,

However it's not at all clear how the current implementation differs to the 2007 reference or where the .phr file is generated from!  Some degree of reverse engineering the C++ source code will be needed!


Metagenomics: Identifying most common sequences

According to the VSEARCH docs, since you have specified --sizeout your abundances have been written into the FASTA headers:

--sizeout
Take into account the abundance annotations present in the input fasta file (search for the pattern ’[>;]size=integer[;]’ in sequence headers). That option is active by default when rereplicating.
Add abundance annotations to the output fasta file (add the pattern ’;size=integer;’ to sequence headers). If --sizein is specified, each unique sequence receives a new abun- dance value corresponding to its total abundance (sum of the abundances of its occur- rences). If --sizein is not specified, input abundances are set to 1, and each unique sequence receives a new abundance value corresponding to its number of occurrences in the input file.


Can exons be located outside of the coding sequence?

Yes, of course. Exons are not limited to the protein coding regions. Many UTRs are in exons. In fact, you even have various cases of UTRs being multiple exons, and being spliced.

What is strange in your file is not so much that you have exons beyond the stop codon, but that you also have them marked as CDS (coding sequence). That isn't possible, no. While there can indeed be exons in the UTRs, those are not coding and shouldn't be called CDS. That's probably just a minor nomenclature issue though and what you are looking at is a spliced UTR.

So, the exons (or the parts of the exons) that fall between the start and stop codons will be translated, while those that fall outside that regions will not be (UTRs).

For example, this is the exonic structure shown for ENST00000617185, one of the transcripts of the human P53 gene:

The boxes are exons and the lines are introns. The colored boxes are protein coding exons, while the empty ones are UTR exons. Note how they are also spliced and that this particular transcript has 3 non-coding exons on the 5' end and one on the 3' end.


Running GSEA or comparable analysis with multiple variable interaction?

I'm not sure I understand "I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables". I would run an association analysis (regression) for module eigengenes and the appropriate interaction term. The (most) significant modules are your candidates. This step is simply regression (in the simplest case a linear model should be sufficient), no GSEA needed. You can then subject the identified modules to GSEA or similar to get an idea of their biological function.


You can use X11 forwarding to run IGV on the headless server where the data is located and simply forward the GUI to your client machine using the X11 protocol. This is a solution I use consistently for my research when working with 20+ tracks that I want to keep on the server.

Download and install IGV for linux from the IGV site on to your server.
Connect to the server using ssh with X11 forwarding enabled

ssh -X user@server
Navigate to the recently downloaded and unzipped IGV folder. Execute the igv.sh script.

./igv.sh

The script launches the IGV session with 4gb of memory. Edit the script to launch with more, 32gb in this edited excerpt.

exec java -showversion --module-path="${prefix}/lib" -Xmx32g \
@"${prefix}/igv.args" \
-Dapple.laf.useScreenMenuBar=true \
-Djava.net.preferIPv4Stack=true \
--module=org.igv/org.broad.igv.ui.Main "$@"


What's the meaning of having different alignments with the same qname in a bam with no duplicates?

I think that you are dealing with two different phenomena:

1) TLEN == 0 reads

From the BAM spec:

  TLEN: signed observed Template LENgth. If all segments are mapped to the same reference sequence, the absolute value of TLEN equals the distance between the mapped end of the template and the mapped start of the template, inclusively (i.e., end − start + 1).14 Note that mapped base is defined to be one that aligns to the reference as described by CIGAR, hence excludes soft-clipped bases. The TLEN field is positive for the leftmost segment of the template, negative for the rightmost, and the sign for any middle segment is undefined. If segments cover the same coordinates then the choice of which is leftmost and rightmost is arbitrary, but the two ends must still have differing signs. It is set as 0 for a single-segment template or when the information is unavailable (e.g., when the first or last segment of a multi-segment template is unmapped or when the two are mapped to different reference sequences).

Note that in each case the mate read appears to be mapped to a different chromosome/contig. So, the mapper can't estimate insert size in these cases, so it just says zero.

2. Multiple reads with same qname

There are various reasons why the same qname can appear multiple times. Again from the BAM spec, one such explanation is that it is a chimeric alignment:

  Chimeric alignment: An alignment of a read that cannot be represented as a linear alignment. A chimeric alignment is represented as a set of linear alignments that do not have large overlaps. Typically, one of the linear alignments in a chimeric alignment is considered the “representative” alignment, and the others are called “supplementary” and are distinguished by the supplementary alignment flag. All the SAM records in a chimeric alignment have the same QNAME and the same values for 0x40 and 0x80 flags (see Section 1.4). The decision regarding which linear alignment is representative is arbitrary.

You may have to do some digging and looking at the different flags to decide what exactly is going on, but my bet is that they are supplementary alignments.


samtools view: writing to standard output failed: Broken pipe

[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting

SAMtools sort has been unable to parse its input, which it thought was SAM (mostly because it couldn't be recognised as another format e.g. BAM).

This is because sed 's/^/LP1-/' is putting LP1- at the front of every line. If the output of samtools fixmate is SAM, then this LP1 is garbling the SAM header lines. If the output of samtools fixmate is BAM (as it probably is), this LP1 is completely garbling the data.

What were you trying to achieve with this sed command?


Create new range based from existing range (python)

Answer to main question:
import pyranges as pr
assert int(pr.__version__.split(".")[2]) >= 93, "pip install pyranges==0.0.93"
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=1, chromsizes={"chr1": 249250621})
gr.Score = np.random.randint(250, size=len(gr))

result = gr.max_disjoint(slack=5000)

# small example that is easy to validate
gr2 = pr.from_string("""Chromosome Start End Score
chr1 1 2 5
chr1 4999 5000 10
chr1 10001 10002 1
chr1 15000 15001 2
chr1 20100 20101 5
""")

gr2.max_disjoint(slack=5000)

# +--------------+-----------+-----------+-----------+
# | Chromosome   |     Start |       End |     Score |
# | (category)   |   (int32) |   (int32) |   (int64) |
# |--------------+-----------+-----------+-----------|
# | chr1         |         1 |         2 |         5 |
# | chr1         |     10001 |     10002 |         1 |
# | chr1         |     20100 |     20101 |         5 |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 3 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Could you please post an example of what data you have and what you want to end up with?

Inserting row after row is going to be slow, whether it is pandas or pyranges. It is best to collect all the data and then create a new dataframe or PyRanges.

pr.itergrs is for iterating over multiple pyranges at the same time. It takes a list of pyranges while you seem to have a collection of tuples. The function iterates over every chromosome/strand pair in the list and for it extracts the data from every pyranges. See the code example at the end.

To iterate over a pyranges you can do for k, df in gr: # ....
Anyways, to do what you want:
import pyranges as pr
import numpy as np

f1 = pr.data.f1()
f1.Start = ((f1.Start + f1.End)/2).astype(np.int32) # int is 64 bit
f1.End = f1.Start + 1

I understand I have not done the addititonal processing you speak of, but if you can post an example I'm sure I'd be able to help you :)

Code example pr.itergrs:
import pyranges as pr
f1 = pr.data.f1()
f2 = pr.data.f2()
for k, dfs in pr.itergrs([f1, f2], keys=True):
    print(k)
     for df in dfs:
         print("df:")
         print(df)

('chr1', '+')
df:
  Chromosome  Start  End       Name  Score Strand
0       chr1      3    6  interval1      0      +
2       chr1      8    9  interval3      0      +
df:
  Chromosome  Start  End Name  Score Strand
0       chr1      1    2    a      0      +
('chr1', '-')
df:
  Chromosome  Start  End       Name  Score Strand
1       chr1      5    7  interval2      0      -
df:
  Chromosome  Start  End Name  Score Strand
1       chr1      6    7    b      0      -


plotting two heatmaps with the same order of genes

ComplexHeatmap is built for plotting side-by-side heat maps with the same clustering - you use the + notation, similar to ggplot2. To use the same example data as @b.nota:
library(ComplexHeatmap)
# First matrix
set.seed(2)
m <- replicate(8, rnorm(26))
rownames(m) <- letters[1:26]

# Second matrix
set.seed(3)
m2 <- replicate(8, rnorm(26))
rownames(m2) <- letters[1:26]

Heatmap(m) + Heatmap(m2)

And the results:

There are lots of examples in the extensive documentation.
EDIT:
Since you have a few genes with high values in the 'T' matrix, they are washing out the colour scheme for the rest of the heatmap. You could set the max to be lower, and accept that those genes will saturate - the resulting heatmap will be less 'honest' about the data, but might capture what you want to see. To change colour scheme requires a call to the colorRamp2 function of the circlize package.
library(circlize)
mycol <- colorRamp2(c(-2,0,2), c("dodgerblue", "black", "yellow"))
Heatmap(m, col=mycol, cluster_columns = FALSE) +
    Heatmap(m2, col=mycol, cluster_columns = FALSE)

This gives:


I would suggest use RALEE—RNALignment Editor in Emacs. It can get for you the consensus secondary structure, you can move left/right sequences and their secondary structures (you can't do it in JalView!), and more.

It's an Emacs mode, so could be a bit hard to start off, but just try, you don't have to use all Emacs features to edit your alignments!

  The RALEE (RNA ALignment Editor in Emacs) tool provides a simple
  environment for RNA multiple sequence alignment editing, including
  structure-specific colour schemes, utilizing helper applications for
  structure prediction and many more conventional editing functions.

Sam Griffiths-Jones Bioinformatics (2005) 21 (2): 257-259.

Fig. You can move left/right sequences and their secondary structures (you can't do it in JalView!)


Paired end sequencing with R1 and R2 of different length. Possible?

As indicated in the comments, yes, you can definitely tell standard Illumina sequencers to sequence mates in a pair to different lengths. This is quite common in single-cell RNA-seq, where R1 is often just a cell barcode and an UMI, so sequencing it to 150 bases would be a huge waste or reagent. One caveat to this is that you can't buy reagents sized for these custom length, just things like 150 base PE or other standards. So unless you end up saving and pooling reagents after such runs you end up paying the same for projects with one shorter read than if you'd just directly sequenced everything to the same length. Given that sequencer runs only take a day or two at the moment, this also doesn't result in that immense of a time savings.


Tool for predicting interactions in the cell

You could try one of these tools to predict protein-protein interactions:

Struct2Net

Given two protein sequences, the structure-based interaction prediction technique threads these two sequences to all the protein complexes in the PDB and then chooses the best potential match. Based on this match, the method generates alignment scores, z-scores, and an interfacial energy for the sequence pair. Logistic regression is then used to evaluate whether a set of scores corresponds to an interaction or not. The algorithm is also extended to find all potential partners given a single protein sequence. Further details about the method are described here.

PRED

Pred_PPI is a web-based system that serves for predicting PPIs from different organisms. This server is freely available to any researcher wishing to use it for non-commercial purposes. Based on auto covariance (AC) and support vector machine (SVM), this tool is capable of predicting PPIs for any target protein pair only using their primary sequences, and assigning an interaction probability to each SVM prediction as well. So the user can use this tool to predict novel PPIs with high confidence.

For protein-RNA interactions, you can try:

catRapid

Through the calculation of secondary structure, hydrogen bonding and van der Waals contributions, catRAPID is able predict protein-RNA interaction propensities with great accuracy (up to 89% on the ncRNA-protein interaction database, NPinter).


How to interpret Mendelian randomization results?

It looks like you have linkage as opposed to linkage disequilibrium, therefore horizontal pleiotropy is an explanation of the calculation.

The calculation appears to be a GLMM around the beta distribution, the residual looks low (which is good). The heterogeneity statistic is significant, you would need to assess whether this is an index of heterozygosity, but it will relate to expectation of Mendellian inhertiance based around the Hardy Weinberg equilibrium.


Control width of samtools tview "snapshot" when redirecting

Partially a shameless plug. See if ASCIIGenome can help you with that.

To review a number of regions that you saved in, say, a bed file, you can either process them in batch or load BAM(s) and bed file and jump from one bed interval to the next using the next -start command.


Find all the bases for given reference position

pileupcolumn.n shows the total number of the reads that cover this basepair (in this case 24793). However, pileupcolumn.pileups iterates over the reads that cover this base pair AND this base pair has a minimum quality in those reads. So it has an additional filter for the reads to be considered. You can modify the quality threshold by setting min_base_quality parameter when you call pileup.
To have a similar filter applied to both coverage and reads, you can use piliecolumn.get_num_aligned() to get the number of segments and pileupcolumn.pileups to get the segments themselves.


What is the per-gene score in BUSCO output?

You are correct, it is the score of the candidate protein against the BUSCO profile hmm assessed by hmmer. If it passes a threshold, it is validated as being that BUSCO gene. The threshold is defined for each profile based on the score of the species used to produce the profile. It can be found in name_odb10/scores_cutoff. And yes, the higher the better.


How to take into account alternative bwa mem mapping when computing coverage

You would want to instead use the -a option with bwa mem to output all alignments. You would be best off filtering these with a custom script to only keep those with a close enough score to the primary hit. samtools depth can't be told to include secondary alignments (it's hard-coded), so you'll need to use something else like bamCoverage or samtools mpileup.


FPKM, FPKM-UQ, TPM or counts: How do I know which kind of unit should I use?

Type of data you need depends on the downstream applications and since you would like to carry out DEA with DESeq2, you would need raw counts (non-normalized).
There are many ways to import/download TCGA data, one such tool, the TCGAbiolinks package gives a nice interface for not only downloading the read count (or pre-processed) data but also associated clinical data.


Impute2 is a great, if dated choice. Another great choice is beagle. The advantage of beagle in this case is that it is not as strict about requiring a recombination map. Without a map, results will be poorer, but perhaps good enough.
Do make sure that your sample size including a reference panel is at least 200 in any case,


You can't run Linux binaries on OSX (or vice versa). Download Mac OSX binaries.


Understanding the "relativeness" in RNA seq experiments

We cannot be sure about this. For all we know, they could have been sequenced to the same depth and there is a biological difference that drives the 1000x (ridiculously large) fold change.
See above

The only thing you can say for sure is that gene B is expressed twice as much as gene A in samples 1 and 2. That's what "relative" means - relative expression of genes within a sample.


How to extract column from a matrix matching the another file (sample file)

Assuming your "master" files looks like this:
df <- data.frame(Name = LETTERS[1:5],
                 Col1 = sample(1:10,5),
                 Col2 = sample(1:10,5),
                 Col11 = sample(1:10,5),
                 Col3 = sample(1:10,5))

  Name Col1 Col2 Col11 Col3
1    A    5    4     3   10
2    B    1    6     5    6
3    C    3    1     9    4
4    D    9    9     8    9
5    E    8   10     1    1

For a sequence of pattern you want to extract, you can use base r function such as grep:
match_pat = c("Col1","Col2")

m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

Col11 Col12  Col2
    2     4     3

And then, extract your columns of interest by doing:
> df[,m]
  Col1 Col11 Col2
1    5     3    4
2    1     5    6
3    3     9    1
4    9     8    9
5    8     1   10

In your example, it will be something like:
match_pat <- c("GTEX.1117F","GTEX.111FC","GTEX.111XS")
m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

master[,m]
```


In SBOL, is there a way to distinguish biological replicates?

Support for distinguishing replicates was added in SBOL 2.2, along with other support for linking data throughout the experiment cycle.  I will explain this in terms of SBOL 3, which is simpler, then give the SBOL 2 translation at the end.
For SBOL 3, all designs, including those for a plasmid, a strain, etc, are specified in terms of the Component class. This is just a "paper" design, however, and may or may not actually exist. When you have a real physical aliquot that you want to talk about, this is represented by an Implementation object: this is your biological replicate.
The Implementation doesn't really have any content of its own; its nature is defined by linking back to the intended design with a prov:wasDerivedFrom link. If you've verified that what you got is what you intended, you can supplement that with a built link as well.
One additional note, however: when you're representing a colony, you're probably not actually representing just a plasmid, but rather are representing a particular strain that has been transformed by a plasmid. The actual data you have may then be one level further removed, since the data has been gathered from a sample under particular growth conditions. SBOL can represent that as well, by describing the sample as a mixture of particular media, inducers, etc., with a transformed strain, which is itself a cell type containing your plasmid.
Here is an example diagram showing just this sort of relationship, of data associated with three replicates of a sample.

You can find this example in the SBOL3 tutorial from COMBINE 2020.
Finally, in SBOL 2, this works exactly the same way, except that instead of a Component your implementation points to a ComponentDefinition or a ModuleDefinition, depending on the nature of the design you're representing. I won't go into that distinction, because it was a bad idea and there's a reason we got rid of it for SBOL 3.


About optimal distance (or dissimilarity) measure in proportion data

I would suggest reading a little on the topic; my understanding is that distance measures can be pretty robust to proportion/compositional data. Here is a relevant textbook chapter:

Other distance measures assume that the data are proportions ranging from zero to one, inclusive.

Among these are indeed Bray-Curtis as well as Jaccard distance (and others). I think that you will want to use one of the more traditional distances rather than a Pearson/Spearman, though they are getting at the same problem. A transformation might help, I suppose, but I would prefer to choose a simple but predictable distance measure rather than a correction on a possibly unsuitable distance that might have unpredictable behavior.
It's a little hard to tell from your question, but it seems that you are wanting to in some way account for indirect correlations using a distance metric. I would not suggest trying to make your distance measure do all this work. Instead, I would use a simple but clear distance method, and then use a multidimensional analysis such as principal coordinate analysis (or similar) to understand the correlation structure of your compositional variables. This is a fairly standard workflow in metagenomics, for instance, which relies heavily on high-dimensional compositional data. I would recommend looking more into metagenomics statistical applications for inspiration if needed.
The choice does also depend to a certain extent on the specific application. In some fields certain methods are familiar and trusted even if there is not a strong empirical/theoretical rule.
You've asked a lot of questions about Pearson/Spearman correlations. Much has been written on the use of correlation analysis on compositional data, most of it cautioning against doing this kind of analysis without special treatment. This leads to difficulties in building meaningful association models for different variables. One method that I'm familiar with in metagenomics is sparse correlation, e.g. SPARCC. But there are newer methods that may be better, for example here and here.


UniProt also provides a BLAST service (https://www.uniprot.org/blast), or you could use the Expasy BLAST too, at https://web.expasy.org/blast.
If you have performed your similarity search at NCBI or elsewhere, you can map the identifiers from that result to UniProt by using the UniProt IDmapping service at https://www.uniprot.org/uploadlists. This can also be used programmatically, cf https://www.uniprot.org/help/api_idmapping


Expected allele frequency distribution of SNVs in real NGS data

I don't have enough experience to answer which probabilistic distribution should be used.

However, this questions also also asks how to estimate parameters of the distributions. If a binomial distribution is chosen, then Heng Li's paper titled "A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data" 1 is probably the definitive one. Section 2.3.1 of that paper describes an EM algorithm for estimating allele frequencies from multiple samples under the assumption of Hardy-Weinberg equilibrium for arbitrary but constant ploidy.


Get gene exons from UCSC browser by HGNC name

You could grab exons via Gencode v27 (hg38) GFF annotations and convert them to BED via BEDOPS convert2bed/gff2bed:

$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "exon"' \
    | convert2bed -i gff - \
    > exons.bed

Then use grep to filter exons against your list of gene names:

$ grep -w -f hgnc.txt exons.bed > exons.filtered.bed


I doubt it, unless you are asking if mapping would be similarly bad for all centromeres. Here are some repetitive structures (probably not centromeric) that I've found in the nanopore reads for "human" sample NA12878, produced by the nanopore-WGS consortium:

These structures are consistent in that they repeat lots of times, but the internal patterns can be quite different. Here are a few more:

Given that centromeres need to be uniquely linkable to a single chromosome, it would make sense to me if the centromere internal structure is unique to each chromosome.

It's possible to have highly-repetitive structures which don't map to each other. While I haven't dug too deeply into the human reads, I've looked at the 5 most-compressible regions in an assembled rodent parasite (Nippostrongylus brasiliensis) genome, and found no internal similarities between them:

One of the issues with assembly from Illumina reads is that these highly-repetitive regions get collapsed into a single repeat (or at best a region up to twice the fragment length). With internal repeat units that have identity of over 98%, assembling the true sequence is very difficult, even with an exact knowledge of paired read separation. Even if this were possible, it can be impossible to correctly place a read because multiple internal units could be identical (or similarly different) to the sequenced read.


install bowtie2 from sources cannot find -ltbb

You need to give the path to the source files for the compiler and the full path to the compiled libraries (all parent directories included) for the linker.  You also need to use either CPATH or CPLUS_INCLUDE_PATH instead of C_INCLUDE_PATH to point to the oneTBB source, since bowtie2 is C++ and C_INCLUDE_PATH will only have an effect for C (see Environment Variables Affecting GCC in the GCC documentation for more info.)
For example like this, inside bowtie2-2.4.1, without making any lasting changes to the shell's setup in this case:
$ ONETBB=$HOME/soft/oneTBB
$ export LIBRARY_PATH=$ONETBB/build/linux_intel64_gcc_cc5.4.0_libc2.23_kernel4.4.0_release
$ export LD_LIBRARY_PATH=$LIBRARY_PATH
$ export CPATH=$ONETBB/include
$ make

After it's built you'll still need the library path to oneTBB defined (via LD_LIBRARY_PATH) so bowtie2 can find it at run-time.
But if you don't have a specific reason to compile from source I'd recommend avoiding it.  You could go with what h.mon suggested and use an existing binary, or use a package manger to install it like zorbax suggested.  Conda has a package for bowtie2 in the bioconda channel, up through and beyond version 2.4.1, and you can use that as a package manager without root access.


Link to UCSC gtexGene track description page for selected gene symbol

You never want to include the hgsid parameter in your url because that is the link to your current browsing session and all its specific contents. So if you send that link to someone else and then they start browsing around that link may not work correctly again in the future if you or the person you sent the link to modify your browsing session.

As to how to create links to hgc pages, the usually follow a somewhat general format:
http://genome-euro.ucsc.edu/cgi-bin/hgc?db=hg19&c=chr2&o=242498135&t=242513546&g=gtexGene&i=BOK

You need the db and g parameters to specify the assembly build and track, and then you need to provide the position in the form of the c,o and t parameters. The last thing you want is the i paramter to specify the item name. The reason the i parameter alone is not enough is because we do not require uniqueness on item names (there may be many copies of an item across the genome for instance).

You can obtain the c, o, t and i parameters by selecting the chrom, chromStart, chromEnd and name fields from the gtexGene table:

$ mysql --user=genome --host=genome-mysql.soe.ucsc.edu -A -P 3306 -e "select chrom, chromStart, chromEnd, name from gtexGene limit 1" hg38
+-------+------------+----------+---------+
| chrom | chromStart | chromEnd | name    |
+-------+------------+----------+---------+
| chr1  |      11868 |    13052 | DDX11L1 |
+-------+------------+----------+---------+


The genes that are missed are probably not official mgi symbols.

You might wanna look them up at mgi:

cd45 -> Ptprc

cd11b -> Itgam

p2ry12 -> P2ry12 (?) (CAPS sensitive?)

glast -> Slc1a3

s100a -> S100a1

pecam -> Pecam1 (?)

I suspect biomaRt is not the problem here, the only better way I can come up with is to download the ENSEMBL annotation (GTF) file, and get the symbol out of this annotation with AWK or in R.

Here an example in R, since you don't want to use AWK:

file <- read.table("Mus_musculus.GRCm38.78.gtf", sep = "\t", comment.char = "#")

exons <- subset(file, V3 == 'exon')

V9 <- file[,"V9"]

ensembl_names <- sub("gene_id ","",V9)

split <- strsplit(ensembl_names, ";")

splitting <- sapply(split, function(x) x[1])

split_symbol <- strsplit(ensembl_names, "gene_name ")

splitting_symbol <- sapply(split_symbol, function(x) x[2])

splitting_symbol2 <- strsplit(splitting_symbol, "; ")

symbol_names <- sapply(splitting_symbol2, function(x) x[1])

ensemble_symbol_big <- cbind(splitting, symbol_names)

ensembl_symbol <- unique(ensemble_symbol_big)

colnames(ensembl_symbol) <- c("ENSEMBL", "Symbol")

write.table(ensembl_symbol, "ensembl_symbol.txt", sep = "\t", row.names = F)


Transcribing DNA exons and append transcripts into a list

I'm afraid none of your code makes much sense. You have written a function that expects an array, and then you are passing it a single letter
for gene in genelist:
    for codon in gene:
        mRNA= transcription(codon)

Then, the function also doesn't make sense. What does next(codon) mean? How would your script know what a codon is or how to move to the next one? And your if codon == 'ATG': is within the if codon != 'ATG': block, so it will never be executed by definition.
I think what you are looking for is something like this, using the same basic approach as your script:
#!/bin/python3
from Bio.Seq import Seq
from Bio.Seq import transcribe

genelist=[]

gene_1= 'A','C','G','G','A','C','T','A','T','T','C'
gene_2= 'G','G','C','C','A','T','G','A','G','T','A','A','C','G','C','A','T','A','G','G','G','C','C','C'
gene_3='G','G','G','C','C','C','A','T','G','A','C','G','T','A','C','T','A','G','G','G','G','C','C','C','A','T','G','C','A','T','T','C','A','T','A','G'

genelist.append(gene_1)
genelist.append(gene_2)
genelist.append(gene_3)

def transcription(ORF):
    mRNA= ''
    foundStart = False
    foundEnd = False
    for i in range(0, len(ORF), 3):
        codon= "".join(ORF[i:i+3])
        if codon == 'ATG' and not foundStart:
            foundStart = True
        if foundStart and not foundEnd:
            cc=transcribe(codon)
            mRNA = mRNA + transcribe(codon)
        if codon == 'TAG':
            foundEnd = True

    return(mRNA)

mRNAs=[]
for gene in genelist:
    mRNA = transcription(gene)
    mRNAs.append(mRNA)
print(mRNAs)

Which returns:
$ foo.py
['', '', 'AUGACGUACUAG']

Of course, it only returns one mRNA since you only have one gene with an ORF in the 1st frame and your script was only looking at the forward strand and the first frame. In reality, ORFs can be in any of the 6 possible frames. You are also assuming that all ORFs will end with a TAG which is perhaps the case in the species you are working in? I don't remember if there is any genetic code that has TAG as the only stop.
If you do actually need to find ORFs in any frame, then I hope the script above should be enough to get you started.


Finding motifs: fasta file with 10,000 sequences

Bio.SeqIO.parse() returns a SeqRecord iterator. To access the sequence of this sequence record, you need to use the .seq attribute so you should update your s with s.seq.


list.of.results <- lapply(seq(INSERT_TOTAL_NUMER_OF_SAMPLES_HERE),function(sample.id){
whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = sample.id,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')[[1]]
})

sig <- do.call(rbind,list.of.results)


Python package or CLI tool outputting mutations of sequence with respect to reference

If it's just a comparison of two sequences that you're looking for, then diffseq from the emboss toolkit should be okay:
https://www.bioinformatics.nl/cgi-bin/emboss/diffseq
I think this tool requires the sequences to be oriented in the same direction. If you discover a lot more variants than you expect, try reverse-complementing one of the sequences.


Finding residues physically interacting in a PDB structure and sequence related to that structure

My reputation isn't high enough to make a comment, but have you considered a modelling software such as PyMol or Chimera? You can execute a command line script with softwares such as these and these have the functionality to predict H-bonds.

In the PDB files there is information regarding the GenBank and UNIPROT info of the protein, and you can find more information here:

http://www.wwpdb.org/documentation/file-format-content/format33/sect3.html.

If you want RefSeq ID, you can use GenBank info you extract and convert by piping into a biomaRt (Bioconducor) script, although maybe there is a more direct way.

As for finding which residues in the sequence are interacting, I'm fairly certain that the software packages like PyMol or Chimera will give you information about the interacting residues.

From what I understand, you will want to filter out only the residues that interact and are on different chains (e.g, Chain A vs. Chain B). This process will require you to work with the syntax of the different software packages when you are scripting. For example, Chimera has a basic primer here:

https://www.cgl.ucsf.edu/chimera/current/docs/ProgrammersGuide/basicPrimer.html

Note that I only gave a reference to Chimera because I am more comfortable with that then other packages and there is likely a more efficient way to do this. With Chimera, you can find the position within the chain of residues participating in the interaction and then you'd store these in a file to do whatever you want later.

I hope this helps.


There seems to be some problem with your "platform specific tempdir" that is used by getGEO's destdir parameter. See ?getGEO for more info on the destdir parameter.
Try
GSE68808 <- getGEO("GSE68808", destdir = getwd())


Unable to read all alignments in stockholm format : UnicodeDecodeError: 'utf-8'

The last line indicates that this is an encoding problem. Try to see if the function you are using to read the data has an "encoding" parameter that you can set to whatever encoding your file is using. If that doesn't work, you can change the encoding of your file by:

with save as using a text editor of your choice and selecting "UTF8" as the encoding. That would be a quick work-around if you have a small number of files.

or

You can use the codecs library to check and convert file encodings within your code. Here is an example.

Warning: I haven't used Biopython in ages but encountered similar problems with Python in the past.


Extract features from fasta sequences and train the classifier

Just a recommendation: Be careful when making your data sets.
P31327_57, P31327_119, P31327_157 are all annotated as "By similarity" to the mouse entry in Swiss-Prot https://www.uniprot.org/uniprot/P31327#ptm_processing, which means there is no direct experimental evidence in the human entry. If you do have publications that report experimental evidence for these PTMs, please don't hesitate to share them with the UniProt team.
I also recommend reading this UniProt help page about using UniProt data in negative datasets if you haven't already done so: https://www.uniprot.org/help/negative_datasets


How to obtain amino acid context at a given position

It turned out that I needed to reinstall OSX to upgrade my R version and biomaRt install.


Assuming that your actual code is this:

for alignment in blast_records.alignments:
    countResults = 0
    for hsp in alignment.hsps:
        if countResults == maxResults:
            countResults = 0
            break
        else:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1

Here you have an outer loop, which is iterating over all the possible database matches, and an inner loop, iterating over the alignments within that database entry. You are currently limiting the inner loop, but not the outer loop, which you probably want to do. Try this:

# Take the first 3 database hits
for alignment in blast_records.alignments[:3]:
    for hsp in alignment.hsps:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)


According to this link, 1=death 0=censor and null=no data.


How to represent trans-spliced genes in GTF?

To my knowledge there's no defined way to deal with that in GTF. GFF3 handles trans-splicing (you'll have to scroll down to "trans-spliced transcript") by giving an individual transcript multiple parents (e.g., ID=some_transspliced_gene;Parent=gene1,gene2). You could use the same methodology with GTF files, but just note that it'll break most downstream programs.


Lambert et al. (2018, Cell. 172(4):650-665) provide a manually curated list of human TFs at this website.


what is the concept of N-channel geneChip?

This is a mere long comment, which may help, but I do not know the details required for the answer.

Several (> 2) samples hybridized in one chip

A traditional microarray has two parts, the mobile target, which is the labelled sample, and fixed to the filter paper, glass or doped silicon is the probe, which is not labelled and is in vast excess. You can label the target fluorescently with one of a variety of fluorophores, or with different radioisotopes, which have different energies. Nobody uses radioisotope microscopy for this except for the most obscure method, say nuclear run off assays, but I mention this just to cover all bases.

Whereas you can buy whatever fluorophore you wish, the detecting machine needs to be set up for it with lasers, beam splitters and filters. Most florescent confocal microscopes use a LED lamp, a splitter and a loads of cubes of filters. In the case of flow cytometry, you can get over 16 fluorescent channels at the same time. However, deconvolution becomes necessary when one channel "bleeds" into an other —I am not sure how one deals with it with a microarray, I suspect before the step to remove regional biases (the step with LOESS regression or similar).

Affymetrix Gene chip (R)

Now, in the case of an Affymetrix Gene chip, I do not know if there is a multichannel machine, but it could be custom build. The spots on an Affymetrix chip are 3-5 µm, so regular microscopy optics still work. Although the silicon base is not transparent like a home made slide or an Agilent microarray (if these are still made). Hopefully, someone else can answer about Affymetrix machinery!

Statistics

I am not familiar with why or how you would analyse a 3+ channel experiment. One guess could be to do a dose-response experiment. But there will be some problems:

The dye flip correction would be an issue for greater than 2 samples
Channel deconvolution may be required for dyes with overlapping spectra.


There's going to be a batch effect between studies, so unless you have some sort of constant group that you can use for between-study normalization you're going to completely muck up your results from using data from multiple studies.

More specifically, this is what packages like sva and the combat() command in it are used to correct. Note that you should go back to the raw .CEL or .fastq files if you want to process multiple studies, in order to at least ensure that everything is processed in the same way.


determining doublets in single-cell RNA-seq

Expected rates of doublets / duplets / multiplets

Fluidigm C1 doublet rate: around 1-5% depending on chip type used.
More information: Fluidigm white paper: Redesign of C1 Medium-Cell 96 IFCs Improves Single-Cell Capture Efficiency (PN 101-3328 B1)
10x Genomics Chromium 3'mRNA-Seq doublet rate: around 1% per 1000 cells captured

Rates in the 'Chromium Single Cell 3’ Reagent Kits v2 User Guide':

Multiplet Rate (%)  # of Cells Loaded   # of Cells Recovered
0.40%   870 500
0.80%   1700    1000
1.60%   3500    2000
2.30%   5300    3000
3.10%   7000    4000
3.90%   8700    5000
4.60%   10500   6000
5.40%   12200   7000
6.10%   14000   8000
6.90%   15700   9000
7.60%   17400   10000

Detection and identification of doublets and multiplets

Experimental

Visual inspection (e.g. Fluidigm C1)
FACS index sorting data can be used to identify doublets.
Mixed species experiment.

Illumina | BioRad White paper Single-Cell RNA Data Analysis Workflow
Zheng et al., 2017: Massively parallel digital transcriptional profiling of single cells

Cell "hashing" with oligo-tagged antibodies (Stoeckius et al., 2017): cells (barcodes) positive for two or more oligos are multiplets.

Computational

Mixed gene expression (e.g. lymphoid and myeloid, or T and B cell)
Ilicic et al. describe an approach for detecting low quality cells (incl. multiplets). They train a support vector machine (SVM) model on a dataset annotated by microscopy inspection. They find, among others, that

"wells containing multiple cells (multiples) show similar expression and noise patterns to broken cells", suggesting "that multiple cells contain a mixture of broken and high quality cells."
difference in proportion of duplicated reads
higher ERCC spike-in ratios

demuxlet (Kang et al., 2017, Nature Biotechnology)
DoubletDetection


What are you expecting these lines to do?

num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate)
...
num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate)

Using the variables that are specified in your code, they are extremely zero-skewed. This is not too surprising given the mutation rate and sequence length.

I expect that increasing the sequence length, or increasing the transition rate, should result in more haplotypes being generated.

After running this through myself, I see that increasing the mutation rate does increase the number of haplotypes, but the skew is retained. This is because each mutation creates a new haplotype, following something like an infinite allele model.

I guess this is not particularly surprising. It's unlikely that two random mutations from a particular haplotype will result in the same haplotype, and therefore any mutations (from any existing haplotype) will generally create something new.


Identification of starting point for the synthesis of a protein

You seem to be interested in identifying potenital coding regions in a DNA sequence. These regions are called open reading frames (ORF) because they are open to translation.
As detailed in the linked Wikipedia article, an ORF usually starts with an ATG and ends with one of the stop codons (TAA, TGA or TAG in DNA space). However, there are many ORFs that are not translated at all, particularly short ones.
Since finding ORFs is an old and established task in bioinformatics there are many algorithms availabe to do this. I would suggest using such an existing tool because they often have additional features and bugfixes. There is a list in the link above.
Estimating the stability and likelihood of a particular amino acid sequence is much more complicated. There are many tools for protein structure prediction that may be worth having a look at.


Remove repeated ALT allele from the REF field of a vcf-like file

I am making a couple of assumptions here:

the delimiter is tab
the reference allele is always present in the ALT column
the reference allele in the ALT column is always in the first position

With those in mind, you can do this using awk as follows:

awk 'BEGIN{FS="\t";OFS="\t"}{$3=substr($3,2)}{print $0}'

Update 1: Using the new VCF style file that @Angel provided, here's how I would do it. This is not significantly different from what @terdon proposed below. The code below will work even if the reference allele is >1 nt.

$ cat test.vcf
#CHROM  POS     ID      REF     ALT
chr4    5708160 .       G       GATAAC
chr4    30280292        .       C       CACGTAT
chr6    166253422       .       C       CCTTTCT
chr11   21931340        .       A       AACACACGC

$ awk 'BEGIN{FS="\t";OFS="\t"}{l=length($4); if ($4==substr($5,1,l)) {$5=substr($5,l+1); print $0} else {print $0}}' test.vcf
#CHROM  POS     ID      REF     ALT
chr4    5708160 .       G       ATAAC
chr4    30280292        .       C       ACGTAT
chr6    166253422       .       C       CTTTCT
chr11   21931340        .       A       ACACACGC


Differential expression between two types of cells of specific marker gene expression

Take a look at the section "Finding Differentialy Expressed Feautres (cluster biomarkers)" of Seurat tutorial.

You have to define your clusters (if not computed automatically) by modifying the ident slot of the Seurat object.

Then you call the function FindMarkers of your Seurat object (to define markers of a given cluster against the rest) or FindAllMarkers (to look for all the combinations)


How to retrieve dataset names from python statsmodels?

Here is the website of the creator of rdatasets I guess with the info you are looking for:

http://vincentarelbundock.github.io/Rdatasets/datasets.html

Here is a short look at it:

Package     Item    Title   Rows    Cols    n_binary    n_character     n_factor    n_logical   n_numeric   CSV     Doc
boot    acme    Monthly Excess Returns  60  3   0   1   0   0   2   CSV     DOC
boot    aids    Delay in AIDS Reporting in England and Wales    570     6   1   0   0   0   6   CSV     DOC
boot    aircondit   Failures of Air-conditioning Equipment  12  1   0   0   0   0   1   CSV     DOC
boot    aircondit7  Failures of Air-conditioning Equipment  24  1   0   0   0   0   1   CSV     DOC
boot    amis    Car Speeding and Warning Signs  8437    4   1   0   0   0   4   CSV     DOC
boot    aml     Remission Times for Acute Myelogenous Leukaemia     23  3   2   0   0   0   3   CSV     DOC
...

I assume the Item column are the identifiers.


downloaded software in Conda environment does not show in ls command

Just type
conda list gzip

OR
conda list | grep gzip

OR if its installed in a given environment
conda activate myenv

or else just activate myenv system depending followed by either of the first two commands.
It is better to use conda-navigate a very nice and under used GUI, and would work in this particular instance. The details are found in an earlier post I made on esentially the same question here together with a pretty picture.
ls will only work if you type ls ~/Anaconda3/bin | grep gzip , but depends on the installation, e.g. you might have installed miniconda.
Having said all that conda is very cool, but whatever you do NEVER use pip or pip3 from this moment onwards EVER, it messes up your conda and you'll need to env your way out of it.


How to get the count of each kmer past 255 using khmer

255 is the default maximum size of a Counttable in khmer. You want to do the following:

import khmer
counts = khmer.Counttable(31, 1e7, 1)  # Feel free to change the parameters
counts.set_use_bigcount(True)

The last line increases the maximum value from 255 to 65535. I don't think there's any way to currently go above that without changing the khmer code (at some point you have to increase the size of the number used in the hash table).


Compute copy number from cases and controls

GISTIC does exactly this, and is very possibly the tool used in your link above.

The input files need to be in a different format than what you currently have, but if you're using a SNP6 array, there are guides out there that will tell you how to get the proper files from your .cel files.


Chunk alignment in a name sorted bam for parallel processing

python will be to slow for this job. Here's a awk solution. One need to sort by read name and take track over the number of reads per chunk. If the number is reached the next chunk can only created of the read name is different.

$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
  BEGIN { part=0; line=n }
  /^@/ {header = header$0"\n"; next;}
  { if( line>=n && $1!=last_read ) {part++; line=1; printf header > part".sam"; print $0 >> part".sam" }
    else { print $0 >> part".sam"; line++; }
    last_read = $1;
  }'

The output will be in sam format.

You can controll the number of reads per chunk by setting the n value for the awk command.

EDIT:

A second version that will output in bam:

$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
BEGIN { part=0; line=n }
/^@/ {header = header$0"\n"; next;}
{ if( line>=n && $1!=last_read ) {print part,line; part++; line=1;}
  print line==1 ? header""$0 : $0 | "samtools view -b -o "part".bam"
  last_read = $1;
  line++;
}


Why do ten rows (Figure_1) correspond to 2 bits (Figure_2) in a sequence logo?

Bits are not frequencies. If a position only contains an A (position 3 for example) then you would need 2 questions (bits) to derive that information without a priori knowledge.

Is it a G or C? if no > then it is a A or T
Is it T? If no then it is an A

Position 1 can be derived by 1 question only:

is it a G or C?  If no then it is an A or T

In this case the number of bits can be calculated using the frequenci)es.

For 4 bases, A only (100%):

4 * -0.25 * log2(0.25) - (-1 * log2(1)) = 2 (bits)

For 8 letters, A only (100%):

8 * -1/8 * log(1/8)- (-1 * log2(1)) = 3 (bits)

Mixed base possibility (A/T 50% each) with 4 bases:

4 * -0.25 * log2(0.25) - (-0.5 * log2(0.5)+-0.5 * log2(0.5)) = 1 (bit)

As mentioned in the comments, you can follow the wiki link for an example.

Or, an example with a little more explanation from the biology SE site.

sequence logo calcuation

Hopes this helps


Reducing a string to a certain amount of the first characters in R

Is there any particular reason why you want to do this via pattern matching in R? This sounds like something that SuperTranscripts and Lace might be useful for:

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1284-1

  Finally, Lace will annotate each superTranscript with blocks either by using the graph structure itself or alternatively using a script called Mobius (packaged with Lace) which infers the annotation from spliced junctions discovered when mapping reads to the superTranscript.


DeSeqDataSet experimental design: column with integer values

Any IDs (e.g. replicate number) don't need to be included in the design formula if it's not a repeated factor in multiple samples.

In this case, you have said that mouse number is related to the number of mice used in the pool (and is a variable factor that is consistent across multiple samples), in which case it should be included in the design matrix. This would also be the case if the same mouse were being used for multiple RNASeq experiments. Converting an integer value to a factor is done by the factor method:

sampleTable$MiceUsed <- factor(sampleTable$MiceUsed);

Other examples of things that might be useful in the design formula would be sex, age, and sequencing run batch (assuming these are not consistent across samples).

With a design formula of something like ~Groups + Mouse, it will also be comparing Groups, it's just hidden in the results. The DESeq2 vignette is really great at explaining what's going on, but as a quick answer you can pull out the Groups results with the contrast function. This will pull out group differences that are independent of mouse number:

results(dds, contrast=c("Groups", "GF", "Wld"))

[first value is the variable name, second value is the numerator condition, third is the denominator condition]


What are negative FSC-A values in flow cytometry data?

This is probably a matter of calibration. If this is not the raw data (actual voltage curves), look at that.

As an area below 0 is not possible, it has to be a calibration issue, either in the machine, or in the software.

the -H plots are about the maximum height of the voltage spike in the side or forward light detector, while the -W plots are about the width (using some definition of flank-detection, so variation prone).
The area is also computed in various ways: using something like (H*W)/c (with c something like 2) or with an actual computation of the area of the spike, but again this has all kinds of predefinitions in it (what's the beginning/end of a  spike, what's the background-level, ...) Some of these predefinitions you can alter in software, some are hardcoded or even part of the functionality of the hardware, so it might be easy to impossible to get rid of the negative values.

The FSC-H values corresponding to the FSC-A values are in your plot, and positive, so the most obvious solution (a bad algorithm that calculated Area as Height times Width times Magic and did not account for negative Heights (which depending on spike definition, you might get)) is out.


Fastq-dump script download X spots or all

You don't need to copy your arguments into variables, you can just use $1 directly. Anyway, this should be enough:

#!/bin/bash
# usage : bash dump_qc.sh <run id> <number of reads>

if [ -z $1 ]; then
  echo "ERROR: at least one argument need to be given!" >&2
  exit 1
fi

if [ ! -z $2 ]; then
    x="-X '$2'"
fi

fastq-dump --split-3 "$x" "$1"


Access workdir defined on command line from within Snakefile

You could query Python for the working directory within the Snakefile, since evidently Snakemake changes the actual working directory of the process to the one specified by --directory.    For example, using a stub Snakefile with just these two lines:

import os
print("Current working directory: " + os.getcwd())

And on the command line:

$ snakemake --directory working-directory
Creating specified working directory working-directory.
Current working directory: /home/jesse/dev/snakemake-wd/working-directory
Building DAG of jobs...
Nothing to be done.
Complete log: /home/jesse/dev/snakemake-wd/working-directory/.snakemake/log/2019-08-02T142009.349338.snakemake.log


GWAS, MWAS, EWAS: what are the (in)dependent variables?

To answer the first part of your question, the dependent and independent variable of X-WAS is kind of arbitrary and dependent on the question you asked. But it gradually becomes a convention in the field after the initial name and concept are accepted by the community. For example, GWAS from the very beginning is written as condition ~ SNP, and there is no doubt about it. Later, people did something like condition ~ gene expression and called it TWAS, which is also widely accepted now. In your case, MWAS is indeed metabolite ~ condition, and the only important thing is to know the convention in the research community so people are on the same page when they talk to each other.

As for why running GWAS using single SNP. In a GWAS, it is typical to use 500K to 7M SNPs. So this is a large number. It is obvious that we cannot include all the SNP in one single model, then the question becomes how many SNPs to include in each model, and how to choose those SNPs. This is very hard to decide. Also, due to link disequilibrium, many SNPs are correlated, which means co-linearity, a thing you want to avoid in linear models. It is not clear how to correctly interpret the effect sizes in a model with many SNPs as well. Therefore, due to feasibility, simplicity, and interpretability, it might be best to perform regression on every single SNP.

Now since we are running GWAS in an SNP-by-SNP fashion, there are a lot of work that needs to be done after GWAS in order to find the true causal SNP. Many fine-mapping methods have been developed to tackle this problem.

Interestingly, there are some multi-variate GWASs methods, but they do not come in the flavor of simple linear models. Instead, sometimes complex Bayesian statistics are used, for instance.


Determine if a gene is mitochondrial or not

You can't determine if a gene is mitochondrial purely from its name. The simplest route would be to use the wormbase biomart and either download the gene names and their associated chromosomes or use the biomaRt package from bioconductor to do the same thing. You can upload a list of the genes you're interested in querying in either case (or just download all of them, and filter).

From your list, only nduo-1, nduo-6, and ndfl-4 are mitochondrial.


Compare alignment quality of multiple sequencing runs aligned against the same reference genome

Qualimap will do this for you.

Go to qualimap.bioinfo.cipf.es
Run qualimap (default params are fine) on each BAM file
Open up the HTML output, and you can read off the %identity (they measure the opposite, i.e. mismatch rate, but 100% - mismatch rate is %identity of course), indel rate, etc.

One thing to watch out for (you don't mention it in your question, but just in case) is that you cannot directly compare Q scores - these are a bit of a mess and calculated very differently in each piece of software.

Unsolicited suggestion: you might also try NGM-LR for mapping MinION data. We've found it beats the others for our data (though we map to a distant reference).


Accessing UCSC genome via ssh results in a validation error

You need to setup ssh keys.  All the commands in the script are going to use ssh to run the commands.  Note the tips here of how to setup ssh keys:

http://genomewiki.ucsc.edu/index.php/Parasol_job_control_system#SSH_keysssh key setup


Negative binomial modeling of RNA-Seq data

The $\mu_{ij}$ in that section of the DESeq manuscript is the expected value of sample $j$ in gene $i$ given its group association (with expected value $\mu_{ij}$). This is why $K_{ij}$ (this is presumably what you meant by $X_{ij}$) is a function of $\mu_{ij}$, rather than $\mu_{ij}$ itself. I think you're just swapping $K_{ij}$ and $\mu_{ij}$ in your head.


Upgrading guppy basecaller on Oxford Nanopore Gridion via PPA

This is an ONT problem, not an ubuntu problem. Keys can expire over time, and need to be updated.
The keys for the Guppy basecaller are not present in the ubuntu repository; you need to get them from ONT. This is indicated by the installation instructions for Guppy (which you've already linked to):
sudo apt-get update
sudo apt-get install wget lsb-release
export PLATFORM=$(lsb_release -cs)
###
wget -O- https://mirror.oxfordnanoportal.com/apt/ont-repo.pub | sudo apt-key add -
###
echo "deb http://mirror.oxfordnanoportal.com/apt ${PLATFORM}-stable non-free" | sudo tee /etc/apt/sources.list.d/nanoporetech.sources.list
sudo apt-get update


How I know which gene is a good predictor in this neural network or not?

Yes, there is a formal solution. You want a heat map of the layers weighting and this is the information sought. It is more difficult to interpret than other supervised learning algorithms, but that is the nature of neural networks.

I know Sci-test Learn (Python), but I don't know R for neural networks, so I can't advise you on the code, only that it will be a few lines of code and a graphical (relative) output. I imagine Tensorflow will have a really brilliant heat map layer weighting.


How to compact variant data to their genes without bias?

Its very easy, just let the ML sort this out for you and that is its advantage, You're thinking of GLM style calculation where you pre-screen the data with bivariate plots, where there needs to be nice Q-Q plots and low residual.
For ML simply include the gene length as one of your parameters along with CpG etc ... and the ML regression analysis SVC, lasso, ridge, random forest will figure the relationship out between gene length and CpG. You do zero, the ML does everything, hence from a statistical point of view purists object because you don't know the relationshiop the ML has deduced between the variables, but you will get regression weights for non-DNN stuff, which will give you some idea of the impact of length.
There is the issue of transformations and that can be complicated, but I'd try untransformed data first. The only disadvantage of this approach is the user will have to input the gene size when they want to check out your training algorithm.


Trying to create a .bam file without the need for a .sam file

This is an error from bwa saying that it can't find your input files:
Here's a checklist:

Check that you have no whitespace after your \ newline breaks. I.e. \ is the very last character
Check that you have no whitespace in any path names
Quote your variables: $REF_PATH -> "$REF_PATH"
Use set -x at the top of your bash script to debug.
In general, read this


What programs account for structural alignment of different parts of distant homologs which have significant structural differences?

TM-align has some options to do this, see for example https://zhanglab.ccmb.med.umich.edu/TM-align/help.html.
TM-align is similar to TM-score but the alignment is purely structural.


How to check what software has been used to generate bam file?

If you're lucky the information is stored in the BAM header in the lines starting with @PG which stores the command lines run on the BAM, use samtools view -H your.bam to display the header.


How to access the patient-derived xenografts (PDXs) repository?

I already see for links when I try to search. Take a look at them below. Also if you can be precise about kind of data are you trying to find, it would be better

https://pdmr.cancer.gov/

http://cdt.northwestern.edu/news/patient-derived-xenograft-repository-now-available-researchers

http://www.epo-berlin.com/epo-tumor-models-xenografts.html

http://www.proxe.org/

http://data-analysis.charite.de/care/

http://www.oasis-genomics.org/

http://www.cbioportal.org/

Hope one of these above will serve the purpose for you.

Update: added some more links that might serve. As for the Cellminer the last update was 2012.


Feature annotation: RefSeq vs Ensembl vs Gencode, what's the difference?

To add to rightskewed answer:
While it is true that:
Gencode is an additive set of annotation (the manual one done by Havana and an automated one done by Ensembl),
the annotation (GTF) files are quite similar for a few exceptions involving the X chromosome and Y par and additional remarks in the Gencode file (see more at FAQ - Gencode).

What are the actual differences between different annotation databases?

They are a few differences, but the main one  for me (and it could be stupid) is
that Refseq is developed by the American NCBI and
the ENSEMBL is mainly developed by the European EMBL-EBI.
Often, labs or people will just start using what is the best known to them (because of a course or workshop) or because they start working with one of the databases with one specific tool and keep with it later.

My lab, for reasons still unknown to me, prefers Ensembl annotations (we're working with transcript/exon expression estimation), while some software ship with RefSeq annotations.

Your lab might be mostly European based people or they might also have read papers like the one from Frankish et al.  Comparison of GENCODE and RefSeq gene annotation and the impact of reference geneset on variant effect prediction. BMC Genomics 2015; 16(Suppl 8):S2 - DOI: 10.1186/1471-2164-16-S8-S2
From the Frankish et al. paper paper:

The GENCODE Comprehensive transcripts contain more exons, have greater genomic coverage and capture many more variants than RefSeq in both genome and exome datasets, while the GENCODE Basic set shows a higher degree of concordance with RefSeq and has fewer unique features.

As for:

Are there significant differences between them today, or are they, for all intents and purposes, interchangeable (e.g., are exon coordinates between RefSeq and Ensembl annotations interchangeable)?

No. I don't think they are great differences between them as that the global picture should stay the same (although you will see different results if you are interested in a small set of genes). However, they are not directly interchangeable. Particularly as there are many versions of Ensembl and Refseq based on different genome annotations (and those won't be interchangeable between themselves either in most cases).
However, you can easily translate most[1] of your Refseq IDs to ENSEMBL IDs and vice-versa with tools as http://www.ensembl.org/biomart/martview for example (there are devoted libraries/API as well like Biocondutor: biomaRt
[1] Most as sometimes, they might be annotated in one of the database but haven't (yet) an equivalent in the other.
EDIT
In fine, even if people tends to keep to what they are used to (and that the annotations are constantly expanded and corrected) depending on the research subject one might be interested in using one database over another:
From Zhao S, Zhang B. A comprehensive evaluation of ensembl, RefSeq, and UCSC annotations in the context of RNA-seq read mapping and gene quantification. BMC Genomics. 2015;16: 97. paper:

When choosing an annotation database, researchers should keep in mind that no database is perfect and some gene annotations might be inaccurate or entirely wrong. [..] Wu et al. [27] suggested that when conducting research that emphasizes reproducible and robust gene expression estimates, a less complex genome annotation, such as RefGene, might be preferred. When conducting more exploratory research, a more complex genome annotation, such as Ensembl, should be chosen.
[..]
[27] Wu P-Y, Phan JH, Wang MD. Assessing the impact of human genome annotation choice on RNA-seq expression estimates. BMC Bioinformatics. 2013;14(Suppl 11):S8. doi: 10.1186/1471-2105-14-S11-S8.


how to filter a multi-VCF for sex-specific genotypes?

Found a great tool with SnpSift!

From the examples: I want to keep samples where the genotype for the first sample is homozygous variant and the genotype for the second sample is reference:

cat variants.vcf | java -jar SnpSift.jar filter "isHom( GEN[0] ) & isVariant( GEN[0] ) & isRef( GEN[1] )" > filtered.vcf


Is it okay to use deeptools bamCompare (SES normalization) for comparisons across different ATAC-Seq datasets?

I would be very hesitant to use SES normalization in a case like this. We recommend that people run plotFingerprint first and see if there is good separation between the samples and only then using SES. I would not expect there to be drastic separation between the samples for ATAC-seq data, so I don't think it will work well. Our standard ATAC-seq normalization is RPGC (aka, "1x normalization"). That should work well for your use-case as well.

As an aside, if you really want to make statistically-backed comparisons between groups I can recommend the CSAW package.