Odds ratio and enrichment of SNPs in gene regions?

The first approach only address the question of how likely are you to end up with the observed over-representation given the MAF distribution.
My suggestion is to use the second approach, but I am not sure if you would call it bootstrap. Bootstrap in general means sampling with replacement to estimate the uncertainty of a parameter, in this case, OR. So even if you resample your full dataset, you get an estimate of how variable your enrichment is, not how likely you are to observe it by chance. so just name it permutation or sampling.
This is what I think might work

Calculate your true OR, qtl snps versus all others for say TFBS
Divide your snps according to MAF into 10 bins / stratas
Within each bin, we permute the label whether the snp is a qtl
Calculate OR of selection / others for TFBS
Repeat 1-3 N times to collect the permutated OR
p value would be number of (permutated OR > observed OR+1)/(N+1)

Below is an example implementation in R, first we have 1e5 signficant snps and they are 5 times as likely to be in TFBS:
snp_is_qtl = rep(c(1,0),c(0.1e6,6.9e6))
snp_is_TFBS = c(sample(0:1,0.1e6,replace=TRUE,prob=c(0.95,0.05)),
sample(0:1,6.9e6,replace=TRUE,prob=c(0.99,0.01)))

We have the MAF for all snps, with those in qtl being different from others and put them in bins.
snp_MAF = c(runif(0.1e6,0.4,0.8),runif(6.9e6,0.1,0.9))

We sort the snps according to MAF to make it easier:
o = order(snp_MAF)
snp_is_qtl = snp_is_qtl[o]
snp_is_TFBS = snp_is_TFBS[o]
snp_MAF = snp_MAF[o]

And we create the bin:
bin_MAF = cut(snp_MAF,10)

Lastly a function to calculate OR:
calOR = function(snps,TBFS){
tab = table(snps,TBFS)
log(tab[1,1])+log(tab[2,2]) -log(tab[1,2]) -log(tab[2,1])
}

So our observed OR (in log) is:
obs_OR = calOR(snp_is_qtl,snp_is_TFBS)
obs_OR
[1] 1.656859

Now we collect the OR under permutation:
N = 100
permutated_OR = sapply(1:100,function(i){
print(i)
perm_qtl = unlist(tapply(snp_is_qtl,bin_MAF,sample))
calOR(perm_qtl,snp_is_TFBS)
})

hist(permutated_OR)

Your observed OR way exceeds the permuted, even after controlling for MAF. Note the computation is a bit heavy because of the large number of snps. So it makes sense to calculate all the ORs (i.e TFBS, promoters..) with each permutation


sort a fasta file containing the Oxford Nanopore Technologies (ONT) header by sequencing start_time ascending

Edit: After doing some reading, if it can be assumed that all the timestamps are ISO8601 and of the same representation of ISO861 (there are a few), then they can be sorted lexicographically, directly, without converting to epoch-time integers. So I have edited this answer to add an answer without conversion steps:

Here's one approach that uses the timestamps directly:

(1) Extract the start times from the headers and sort them to another file, in ascending order:

$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa | sort > starts_iso8601.sorted.txt

(2) Use that file of ordered ISO8601-formatted starts to write out FASTA records in your sort order:

$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

Quite a bit cleaner and faster.

You might separate the timestamp extraction and sort steps into two commands and inspect the intermediate result, just to ensure that you have valid strings. Or use the approach below to convert the timestamps, which will catch bad input.

Here's a second way that may add unnecessary work. As an advantage, however, it adds some error checking by converting the ISO8601-formatted timestamp.

If the string you pull out of the header is not a valid timestamp, the first conversion step (ISO8601-to-epoch) will fail and you will then be able to catch invalid input pretty quickly, which can add confidence that your output is correct.

(1) Extract the start times from the headers:

$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa > starts_iso8601.txt

(2) Convert the ISO8601-formatted times to time-since-epoch times:

$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"$line" +%s; done < starts_iso8601.txt > starts_epoch.txt

(3) Sort the epoch times — which are just integers — to another file, in ascending order:

$ sort -n starts_epoch.txt > starts_epoch.sorted.txt

(4) Convert sorted epoch starts back to ISO8601 format:

$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"@$line" -u +"%Y-%m-%dT%H:%M:%SZ"; done < starts_epoch.sorted.txt > starts_iso8601.sorted.txt

(5) Use that file of ordered starts to write out FASTA records in your sort order:

$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

Note: Use gawk, gdate, gsort in place of awk, date and sort if you're using OS X, Homebrew, and GNU core utilities.

I'm showing the steps separately to explain the general procedure.

Many or all of these steps can be glued together in a Unix pipeline to save the time, memory and disk space costs of writing intermediate files.

Also, I'm not doing much error checking here. You might look at cases where headers do not contain start_time, or times are not all in ISO8601, or other cases not anticipated here.


It means the distribution is heavily skewed and has a 'long-tail', thus the variance is much greater than the mean. In other words most values are around zero, with a small number of values having extreme values.
[
You need to check the goodness of fit because that is saying the deviation from a power-laws distribution is just significant. So it ain't quite power-laws, but the general gist of a 'long-tailed distribution' is correct.
Its a good result and a nice analysis, but you are not quite there at present (needs to be >0.05). Check the negative binomial distribution.

I've seen you data, it looks very good and the power series is not a bad explanation of the data. There are two possibilities for the goodness of fit being 0.035

The goodness of fit is low because of 'outliers' at the beginning of the distribution
An alternative distribution, such as a negative binomial distribution fits better
Check it is 0.035 not 0.35 because that changes everything.

The best way to do this is to plot the power series expected distribution onto the graph and look where the fit 'fails'.
If it is not fitting at the beginning of the distribution thats is about outliers and easily solvable. If it is not really capturing the shape then its probably a slightly different distribution.
If its point 1 then just remove the very first point in the graph and try again and check the p-value. Or take the mean of the clump of three points next to the first point and try again (or just remove those three point). Don't do any more than that and the technical word is "outlier" analysis. In your write-up you must explain the issue around any 'outliers' and that is perfectly acceptable. If it is a power-series distrupted by an outlier(s) there is formal analysis to demonstrate an 'outlier' using outlier analysis if any one ever asks you.
If its point 2. then use e.g. negative bionomial or the gamma distribution (alpha <1) and check the goodness of fit.
You want a goodness of fit > 0.05 and preferably >0.1 to be sure its the right distribution. There are formal statistics to do this, but that is the main one.
My guess is there are outliers at the beginning because your sample size is not large enough to identify the variation correctly, hence the fit is low. In 'long-tails' its easy to have a 'low sample size' because of the way they work. Thus in this case mapping a distribution is a perfect way to compensate for the sampling.


How to obtain .bed file with coordinates of all genes

Via Gencode and BEDOPS convert2bed:

$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

You can modify the awk statement to get exons, by replacing gene with exon.

BEDOPS: https://github.com/bedops/bedops

This is based off an answer I wrote on Biostars, which includes a Perl script for generating a BED file of introns from gene and exon annotations: https://www.biostars.org/p/124515/#124522


Glocal\semi-local\Hybride Globale-Local alignment with Python

As suggested in the comments, the built-in BioPython pairwise aligner should work for this. You can find a short sequence match in a longer sequence using a global alignment with the right score parameters (query_end_... = 0). The other scoring parameters are values that I've found to be good for aligning similar sequences.
Example code:
from Bio import Align

aligner = Align.PairwiseAligner()
aligner = Align.PairwiseAligner()
aligner.match_score = 5
aligner.mismatch_score = -9
aligner.mode = 'global'
aligner.target_internal_open_gap_score = -12
aligner.target_internal_extend_gap_score = -3
aligner.target_end_open_gap_score = -12
aligner.target_end_extend_gap_score = -3
aligner.query_internal_open_gap_score = -12
aligner.query_internal_extend_gap_score = -3

# Set the query_end penalties to zero to effectively implement a 'glocal' alignment as described in the question
# These are the default values, so you could skip the next two lines,
#    but this is relevant to the question, so I'm including them in the example.
aligner.query_end_open_gap_score = 0
aligner.query_end_extend_gap_score = 0

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATCCG")

for a in alignments:
    print(a)
print()
print()

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATACG")

for a in alignments:
    print(a)
print()
print()

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATACCG")

for a in alignments:
    print(a)


Yes, you end up placing filters where things become inconsistent or similar to background. In the example, which appears to be from the Seurat tutorial, there's a downward plunge in nGenes vs. nUMI at around 200 genes and few believable cells with more than 2000-2500 genes, so the tutorials sets filters around there. The 0.05 mitochondrial contamination is a more arbitrary round number and is less supported by the plots (one would use more like 0.06 or 0.07 if one based this solely on the plots).


Another option would be to adjust the quality of the primer bases to be below the threshold that your variant caller uses.  In my case I have "masked" the primer bases by dropping the quality to 0.  This prevents the variants from being called based that are synthetic but if you have an overlapping amplicon it will call bases from those reads.


Why did expression based subtypng of breast cancer gain much more acceptance than others

In two words: incidence and funding
I'm not an expert on this topic, but I assume it has something to do with the incidence of breast cancer itself:

Breast cancer is the most common cancer in American women, except for skin cancers. Currently, the average risk of a woman in the United States developing breast cancer sometime in her life is about 13%. This means there is a 1 in 8 chance she will develop breast cancer. This also means there is a 7 in 8 chance she will never have the disease.

This incidence most likely amplifies how big of a problem society/funders/governments consider breast cancer to be. For example, a study performed by Northwestern University’s Feinberg School of Medicine and Robert H. Lurie Comprehensive Cancer Center showed that breast cancer research is by and large one of the most well-funded areas of cancer research.

This higher amount of funding leads to more and better resources and research into every aspect of breast cancer, from identification to treatment. I would wager better funding has had the effect of increasing the level of research performed in the usage of sub-typing for identifying people at risk of getting breast cancer. Given positive results and lots of literature backing up such a method as useful, its levels of acceptance in clinical settings increased.


Demultiplexing FASTQ file without index information

I would suggest looking at sample metadata:

Library:

Name: TruSeq RNA

Instrument: Illumina HiSeq 2000

Strategy: RNA-Seq

Source: TRANSCRIPTOMIC

Selection: PolyA

Layout: PAIRED

You can see that it's paired. You can look at the spot descriptor at the page you link to see more information about the design.
There is no multiplexing information available there. I would instead look at the overall project and its metadata to find the different replicates. It looks like there are 2 samples x 3 replicates each, deposited independently. So it is likely to already be demultiplexed, and you can just download them individually.
"Unstranded" isn't really a meaningful term for this, if I understand correctly. Forward and reverse read pairs don't correspond to forward and reverse on the chromosome, but to the way that sequencing-by-synthesis happens on the machine. So the mapping orientation is not really relevant.
I infer that you are working from a single FASTQ file. If you are working from a single FASTQ file for paired-end out of SRA, then the FASTQ file is not yet resolved into F+R reads. I would suggest rerunning the FASTQ conversion with --split-3 or similar option.


How can longest isoforms (per gene) be extracted from a FASTA file?

While the solution from https://bioinformatics.stackexchange.com/users/96/daniel-standage should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:

#!/usr/bin/env python
from Bio import SeqIO
import sys

lastGene = None
longest = (None, None)
for rec in SeqIO.parse(sys.argv[1], "fasta"):
    gene = ".".join(rec.id.split(".")[:-1])
    l = len(rec)
    if lastGene is not None:
        if gene == lastGene:
            if longest[0] < l:
                longest = (l, rec)
        else:
            lastGene = gene
            SeqIO.write(longest[1], sys.stdout, "fasta")
            longest = (l, rec)
    else:
        lastGene = gene
        longest = (l, rec)
SeqIO.write(longest[1], sys.stdout, "fasta")

If you saved this as filter.py, then filter.py original.fa > subset.fa would be the command to use.


Experimental Design for Differential expreression analysis

You're really running out of degrees of freedom and the only actual replicate is a single sample from another patient, so please take the results with a huge grain of salt.

You can use a design of ~cellType * coculture where cellType has levels NOF and CAF and coculture has levels yes and no (ideally you'd block by patient, but that's not an option for you).


How to retrieve the yeast locus tag from uniprot via sparql?

Ok, figured it out by looking at the scheme uniprot provides:

PREFIX up:<http://purl.uniprot.org/core/>

SELECT ?protein ?yeast_id
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein up:encodedBy ?gene .
    ?gene up:locusName ?yeast_id
}

So one starts at the accession number in the middle and then uses the path to the actual locusName.

Also the path via the database works fine, I just did not use the seeAlso correctly (protein was only bound to db but not to db2, that's why it returned all entries):

PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT ?protein ?db ?db2
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?protein rdfs:seeAlso ?db2 .
    ?db up:database <http://purl.uniprot.org/database/SGD> .
    ?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

}

This will return the desired output (one would still have to process the db2 entry, of course):

protein db  db2
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446  http://purl.uniprot.org/eupathdb/FungiDB:YOL086C


Can a customized GRCh38 .gtf file be used with any of the GRCh38 released patches?

You can see the exact assembly by looking in the README file in Ensembl's genome download page. ftp://ftp.ensembl.org/pub/release-89/fasta/homo_sapiens/dna/README

As you can see the current assembly is GCA_000001405.25. Ensembl 79 uses version GCA_000001405.17. You can expect some differences between assembly versions.

If you take a look at https://www.ncbi.nlm.nih.gov/assembly?term=GRCh38&cmd=DetailsSearch

You can see that 20,466,394 bp have been added since 17 -> 25. This is only about 0.6% change so I would guess there is not that much difference between these versions. However, if I were you I would use the corresponding version of Ensembl. Even better would be using the latest version of Ensembl, as not only the genomes primary sequence but also the positions of genes can change between versions (https://doi.org/10.1093/bib/bbw017).


Is BWT based aligner suitable for any types of alignment task?

BWT based algorithm is just a component of a full-pledge aligner. To that end, it is applicable to most alignment tasks. However, for long-read alignment and metagenomes, we don't often use BWT based algorithms as other methods are faster.


Parse BAM by insertion size and get genomic coordinates

You basically have the right framework already in your pysam code, which I took the liberty of tidying up a bit. You just need to add a bit:

for read in bamfile:
    if not read.is_unmapped:
        cigarLine = read.cigar
        interesting = False
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1:
                if cigarLength > 100:  ## Arbitrary value
                    interesting = True
                    break
        if interesting:
            print("{}:{}-{}".format(read.reference_name, read.reference_start, read.reference_end))

Unless you have long reads, this will be close enough for IGV.


Multithread fastq processing with kseq.h in C++11?

You might take a look at my FQFeeder library (https://github.com/rob-p/FQFeeder), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the parsed reads.  The library contains an example program to demonstrate how to use it with multiple threads.  It's worth noting that kseq is already very fast, so multi-threading only really makes sense if worker threads are doing non-trivial work on each read so that you can then make the parsing and processing concurrent (have the processing keep pace with kseq).

I'll note that my library relies on a header-only MPMC implementation and so technically isn't "only using the standard library", but you're already relying on kseq.h, so omitting third-party libraries must not be a hard constraint.


differential gene expression complex design no replicates

Although it is not recommended to use no replicates, in the edgeR manual they give some advice on how to go on with no replicates design. See page 21 of user guide. You can e.g. estimate a BCV value. Of course it is still a trick and not sound statistics.


bed file with N regions of GRCh38 reference?

# if you have seqtk installed, skip the following two lines
git clone https://github.com/lh3/seqtk
cd seqtk && make
# the main command line
./seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed

Option -n sets the min stretch length. Just use -p as is. It is a bit complicated to explain what it is doing.


but I would say the pausing index is something that is invented, commonly used, but with no clear indication of what is the expected distribution or how it can be subjected to test like you've mentioned.

I assume you would like to get your work published, so going by the papers with IF > 5 such as this, and this, and this, the clear strategy is use a boxplot to show the difference between conditions, and apply a statistical test such as mann whitney u or t.test to say, there is a difference.

Let's try this, sorry not very good in R,using your dataframe which i named x:

library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)

df = x %>% select(contains("pausing") | contains("Transcript")) %>%
rename_at(vars(1),substr,1,10) %>%
pivot_longer(-Transcript,names_to = "condition") %>%
mutate(condition=case_when(
grepl("Control",condition) ~ "Control",
grepl("PTPN6",condition) ~ "PTPN6_KO"
)) %>%
group_by(Transcript,condition) %>% summarize_all(mean)

In the data frame above, I have averaged the pausing index between replicates (not very good but then again, most papers have only 1 replicate), and group them by condition.

# A tibble: 12 x 3
# Groups:   Transcript [6]
   Transcript   condition value
   <chr>        <chr>     <dbl>
 1 NM_001011666 Control   2.84
 2 NM_001011666 PTPN6_KO  1.86
 3 NM_001301039 Control   1.48
 4 NM_001301039 PTPN6_KO  1.14
 5 NM_022114    Control   1.03
 6 NM_022114    PTPN6_KO  1.31
 7 NM_173803    Control   0.834
 8 NM_173803    PTPN6_KO  1.37
 9 NR_132775    Control   1.08
10 NR_132775    PTPN6_KO  1.33
11 NR_160728    Control   0.818
12 NR_160728    PTPN6_KO  0.892

Then we plot:

p <- ggboxplot(df, x = "condition", y = "value",
                color = "condition", palette =c("#c72e29","#016392","#be9c2e","#098154"),
                add = "jitter")
p+stat_compare_means(comparisons=list( c("Control", "PTPN6_KO")))


Getting Unique Identifier List for GEO Datasets NCBI

You can use Entrez Direct for this. The following returns Unique Identifiers which are just bare integers.
$ geo_query='"Expression profiling by high throughput sequencing"[DataSet Type] AND ("Homo sapiens"[Organism] OR "Mus musculus"[Organism] OR "rattus norvegicus"[Organism]) AND ("2020/01/01"[PDAT] : "3000"[PDAT])'
$ esearch -db gds -query "$geo_query" | efetch -format uid > gds_results.txt
$ wc -l gds_results.txt
9981 gds_results.txt
$ head -n2 gds_results.txt
200134092
200120931

Instead, if you are looking for a way to get the GSE accessions, you can use the built-in xtract command to parse the XML returned by esummary as follows:
$ esearch -db gds -query "$geo_query" | esummary | xtract -pattern DocumentSummary -first Accession > gse_accs.txt > gse_accs.txt
$ wc -l gse_accs.txt
9981 gse_accs.txt
$ head -n2 gse_accs.txt
GSE165829
GSE165824


Extract LINE-1 methylation status from WGBS data

After a more extensive search, I found a useful answer that I'll share in case anyone else wonders how to do this. I'm using MethylKit package for differential methylation expression and it is quite easy to integrate with Genomic Ranges and rtracklayer in R to extract particular element methylation differences:

element.of.interest <- import.bed("FILENAME.bed") #download custom track from Genome Browser
element.meth <- diffmeth[as(diffmeth, "GRanges") %over% element.of.interest]

where "diffmeth" is the resulting MethylKit object from the differential methylation analysis.


Generate an approapite output from DNaSP6 to Arlequin program

Yes, your question makes alot of sense and its a good question because hits at the heart of differences between population genetic analysis and phylogenetic analysis. However, when you speak about gaps you are referring to 'indels' (insertion/deletion).

Indels The conservative approach is to simply delete sites with gaps. Arlequin uses a lot of maximum likelihood solutions and gaps would automatically be deleted in an Arlequin max. likelihood analysis, because it is the output per alignment per matrix. If the matrix doesn't have a 'gap'.

Indels and pop gen Thats tricky, it is better to remove them. HOWEVER the indel may contain underlying population signal, if its a maximum likelihood analysis (pop gen does this too) I suspect they will be ignored. If its pairwise, such Fst including indels might give you a better answer because there is more polymorphisms to assess population structure. So if you're doing a pairwise analysis, I'd do it both ways and then report the methodology correctly if there are differences in results and simply explain the theory.

Invariable sites This is different from an indel. An invariable site is a site which carries no phylogenetic signal, it carries information but a phylogenetic analysis could draw any conclusions because all the nucleotides/amino acids are identical at a give alignment position.

The output you obtain in allele frequency will be different and here's the rub ...

Phylogenetics analysis Arlequin is really pop gen, but it will contain some phylogenetics programs. If you see anything with 'Jukes Cantor correction', 'Kimura 2 (or 3) parameter model', 'HGW model', '"REV" model' (eek I can't remember its name any its a reverible 6-parameter model of nucleotides),  you are dealing with phylogenetics. Keeping invariable sites is normal here because the rate matrix is calculation "per alignment site".

Population genetics analysis Fst, Fis, linkage disequilibrium statistics, REMOVE THEM. Why? Because you want SNPs (single nucleotide polymorphisms). The invariable sites will get in the way of the statistics. For example is the locus recombines and you want to assess a linkage disequilibrium statistic (recombination), invariable sites will recombine, but the algorithm can't see this because the nucleotides are all the same. The only thing the algorithm can see and analyse are SNPs, which 'representative' of all the surrounding invariant sites. Thus removing invariants is essential. Arelquin might do this automatically (but I don't know the preparatory algorithms), but if you submit multiple loci Arlquin will want to know which are different loci. It will then recode the loci automatically if a single locus comprises more than one SNP. We are going into the underlying theory of linkage disequilibrium here and I'm not sure if say Fst/Fis is particularly bothered with this or is just concerned with the populations of genetic diversity.

Summary So gaps ... if they are genuine indels and not just 'N', ie indeterminant nucleotide calls, remove them for phylogenetics. If its pop gen do the calculation both ways and assess the output, if there is no difference report the 'gap/indel removal' result.

Invariant sites .. for pop gen you remove them but Arlequin might do this automatically to prepare a data set for analysis (I don't know). For phylogenetics you never remove invariant sites.

Goodluck


My understanding is that PBJelly is the tool intended to solve this problem.
This person seems to have used it starting with Illumina scaffolds. I have not personally used it but have heard that it works pretty well.
Edit: note that the PacBio ecosystem is not super easy to use from scratch. However, they have made pretty good efforts to facilitate setup using conda, as for instance documented here. This person shows a conda recipe for pbjelly.


How to specifically select cluster on Featurescatter on Seurat?

You can use the cells argument to only display the data for a specific subset of cells.

For example to extract the values for a cluster 2:

# The cluster identity is stored in the meta.data
# In this example there is a clustering column called "RNA_snn_res.0.7"
# but any metadata column can be used

clust2.cells <- rownames(Seurat.data@meta.data[Seurat.data$RNA_snn_res.0.7 == 2, ])

FeatureScatter(Seurat.data, feature1 = "gene_A", feature2 = "gene_B", cells = clust2.cells)


How can I download the SNP genotyping information from GTEX?

There are two types of data you can get from GTEx. The first is summary statistics, basically the slopes and p-values you see in the boxplot. You can also get expression levels for all genes measured in GTEx.

However, in order to plot boxplot by yourself, you also need genotype data. This is not publicly available because of privacy issues. In order to get it, you have to apply to it on dbGaP. I think you can find relevant information on the GTEx website.


Different outputs of computeMatrix depending on binsize

Please don't send things like this to the deepTools mailing list and post them here.

As I mentioned on the list, what you're likely being surprised by is that computeMatrix is not summing values, but taking the average in each bin. You can change this with the --averageTypeBins option (the default is mean, you presumably want sum).

Update: For those wondering the cause of the confusion was due to not using --missingDataAsZero. In other words, including or not including NAs will lead to different results, since different numbers of values will be included.


The NCBI BLAST can be used with the PDB DB (which NCBI has). The PDB codes are stored as 4 letter codes underscore chain, e.g. 1GFL_B. The catch is segment identifiers —but generally they are the same peptide so shouldn't be an issue.
You can search specifically the PDB DB in NCBI (not the RCSB PDB) by setting the database to PDB.


What are the state-of-the-art cell-type RNA-Seq deconvolution methods?

I see a lot of people using xCell. There should also be papers systematically compare many different methods, which might interest you.


You could use BEDOPS bedmap to map reads to introns, using 1) the --count operator to do counting of reads overlapping by your criteria; and, 2) the --indicator operator to do a true/false operation, where reads are contained entirely within the intron.

For instance, to count reads that overlap introns by at least 25 bases, use --count to count and --bp-over 25 to specify the minimum overlap threshold:

$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

To filter for introns where reads are contained entirely within the intron, use --indicator to get introns that meet the criteria of full containment of reads, via --fraction-map 1:

$ bedmap --indicator --echo --fraction-map 1 --delim '\t' introns.bed reads.bed | awk '($1==1)' | cut -f2- > answer.indicator.bed

Once you have answer.counts.bed and answer.indicator.bed, you can compare those two results with bedmap --indicator --exact to see where there are commonalities (true/1) and differences (false/0), and increment the final count where there are identical introns, e.g.:

$ bedmap --indicator --echo --exact --delim '\t' answer.counts.bed answer.indicator.bed | awk '{ $NF += $1; print $0; }' | cut -f2- > answer.final.bed

If your reads are in a different format, i.e. BAM, you can use bam2bed in a process substitution, replacing reads.bed with <(bam2bed < reads.bam), e.g.:

$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed <(bam2bed < reads.bam) > answer.counts.bed

Etc.

You can also just convert separately, so that you only need to do it once:

$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

Etc.


How can I find the relevant pathway map from gene-gene or protein-protein interaction list?

See a project like Plant reactome, where they infer the pathways of all sequenced plant genomes using orthology to well annotated species.
http://plantreactome.gramene.org/


You can do this with the reutils package, which provides an API to NCBI's E-utilities. Here's an example for your specific question:
install.packages("reutils")
library(reutils)

# Get universal identifier
uid <- esearch("NP_000029", db = "gene")
# Fetch summary
sm <- esummary(uid, db = "gene")
# Extract specific tag
sm$xmlValue("//Name")


Do I need to fill in missing side chains and loops before protein-ligand docking for VLS?

In general
Most if not all docking applications will tolerate a jump.
Furthermore, you really don't want to dock anything to the flexible regions even if you are using a program that does properly protein flexibility.
Furthermore, if your protein adopts a closed conformation upon native substrate binding (not really the case here), you should start with that for better docking results.
Cytochrome P450
So most Cytochrome P450s are strange though and in many those loops play a large part in their catalysis. I am on a paper under review (bioarxiv) that deals exactly with this in a different cytochrome (P450-BM3, the biotech one). Namely even if the active site is fairly large for an active site in the core of an enzyme, the ligand (testosterone) did not really rotate around and instead the specificity of oxidation was controlled by how the ligand got in. So those crystallographically unresolved loops act as gatekeepers for what can bind. In terms of VLS, one aims for blocking off the cofactor, i.e. the heme. So you can do docking without worrying about them in a first instance, but there is a possibility that a lot of active drugs instead interact with the loops.
Certainly, dynamic undocking and other MD methods may be required a few iterations in, to see why some candidates failed and why others worked.
An interesting unofficial further possible spin is the fact that solvent-exposed sequential residues Asp-Pro can hydrolyse under acidic conditions, in a mechanism similar to isoaspartate formation, and several cytochromes have these on the entrance loops and have weird kinetics and smeary SDS-PAGE gels.
Lastly, you are totally correct in your observation. The water binding the heme in the crystal is very special. It is taking the place of oxygen. That could be modelled as is, changed to a pre-encounter O2, or event changing the heme to a heme bound to a radical oxygen —note that some forcefields (CHARMM, Rosetta Ref2015) do not like this weird valency iron.


Removing stop codons from alignment file for dn/ds calculation

I would align your sequences with the software available here: https://github.com/veg/hyphy-analyses/tree/master/codon-msa


Read and write FASTA files with more information than id and sequence

The description field in the SeqRecord object has the information you are looking for:

>> from Bio import SeqIO
>> s = SeqIO.read('genome.fasta', 'fasta') # single sequence fasta file
>> s.id
k119_5
>> s.description
k119_5 flag=0 multi=141.0706 len=473

Edit: as an aside, if you write a SeqRecord object using SeqIO's write method, the FASTA header should contain the whole original header.


Optimal design for an RNA-seq experiment with ERCC RNA Spike-In Control Mixes

I have 6 conditions and 3 replicates in each condition.

If including more biological replicates has a low experimental cost, you should be using 6 replicates per condition (and fewer reads per condition, i.e. maintaining the total number of reads) in order to pick up the majority of differentially-expressed genes. For that, I'd recommend using 3 lanes of a flow cell, with 12 dual-indexed barcodes (with the same index on each end) for each lane. Doing dual indexing will protect a bit against index switching and make the results more robust.

  randomly assigning mixes 1 and 2 among samples. The latter is not very clear to me. Should I use the same mix for all replicates within a condition or should I randomise really among all the samples?

Randomly assign mixes among and within the conditions. Make sure that the mixes are equally represented (as much as possible) across the conditions, but that each condition has at least one of each mix. This can be done using something like the following in R:

numReplicates <- 7;
numEachMix <- ceiling(numReplicates / 2);
assignments <- sample(rep(c("mix1","mix2"),numEachMix))[1:numReplicates]

For working out the level of fold change this is reliably replicated (or "statistically significant"), you can split the samples into mix 1 and mix 2, then apply the same algorithms for ERCC differential expression as with the non-ERCC genes.

However, as ERCC is typically added in proportion to the amount of RNA (or, if you're lucky and have cell-sorted populations, possibly the number of cells), all it's really testing is the accuracy of sample prep.


ipfam or pfam domain-domain interactions

I think you're looking for ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/database_files/pfamA_interactions.txt.gz. Based on the file's name, that's a list of domain interactions. The file has two columns, so presumably each line represents an interaction.

Here's what I get for the Death domain (PF00531):

$ zgrep PF00531 pfamA_interactions.txt.gz
PF00036 PF00531
PF00531 PF14786
PF00531 PF00531
PF00791 PF00531
PF01335 PF00531
PF13499 PF00531
PF13833 PF00531
PF14786 PF00531

Those are 8 interactions, which matches what the PFam link you gave in your question shows. If you want to convert the accessions to something more informative, you could try Pfam's REST API and do something like:

$ zgrep PF00531 pfamA_interactions.txt.gz |
    while read a b; do
    name1=$(wget -O - "https://pfam.xfam.org/family/${a}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    name2=$(wget -O - "https://pfam.xfam.org/family/${b}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    printf "%s\t%s\t%-20s\t%-25s\n" "$a" "$b" "$name1" "$name2";
done
PF00036 PF00531 EF hand                 Death domain
PF00531 PF14786 Death domain            Tube Death domain
PF00531 PF00531 Death domain            Death domain
PF00791 PF00531 ZU5 domain              Death domain
PF01335 PF00531 Death effector domain   Death domain
PF13499 PF00531 EF-hand domain pair     Death domain
PF13833 PF00531 EF-hand domain pair     Death domain
PF14786 PF00531 Tube Death domain       Death domain

Obviously, if you need anything more fancy than that, it's time to write a proper parser using Perl or Python or any other language that can understand XML.


Set 1

read 1 6523432: A - 6710432 T
read 2 6523432: C - 6710432 G
read 3 6523432: A - 6710432 T
read 4 6523432: C - 6710432 G

I understand this as showing the location of 2 variants, supported by 2 reads each. Now, if the reference has an A at position 6523432, your sample could be heterozygous for an A->C variant at position 6523432. If the reference has a C at that position, then your sample could be heterozygous for C->A. If the reference has something else again, for example G, then your sample could be heterozygous for G->A/C. Similarly for location 6710432.
Set 2

read 1 6523432: A - 6710432 T
read 2 6523432: A - 6710432 T
read 3 6523432: A - 6710432 T
read 4 6523432: A - 6710432 T

In this case, if the reference has an A at position 6523432, the reads do not support a variant. However, if the reference has anything else, then the reads would support a homozygous variant e.g. G->A. And the same for position 6710432.

Now, in all of these cases, you cannot know whether a variant is real without taking into account the number of reads that support it, the quality of the alignment, and the quality of the reads. That's why we do this using dedicated variant callers, software that is designed to take all these factors into account and decide whether to call a variant at that position.

But in all cases, the variant must be defined with respect to a reference genome. You cannot call a variant just by looking at reads without any knowledge of a reference. A "variant" just means "not what is found in the reference". It is largely a matter of random chance whether something is a "variant" or "normal". If the individual(s) whose genome was taken as a reference has the same sequence, then you call it normal. If not, then you call it a variant. There are many cases where the standard hg19 or hg38 genomes have a "normal" residue which is actually far less common than the known "variant". There's nothing special about the reference genome, after all, it's just the one we chose to call a reference.


The best solution I could find was to use the \shaderegion function on all the sequences so it ended up looking something like this:

\shaderegion{1}{CG}{Red}{OrangeRed}

\shaderegion{2}{CG}{Red}{OrangeRed}

...

Then I cast the rest of shading similarity to grey so the CG pairs were more easily distinguishable.

\shadingcolors{grays}

Two examples one with \fingerprints and one regular.


Which correlation method to compute the correlation score between different clusters of Sc-RNAseq data?

Disclaimer: This is not a complete answer but more a comment, posting as an answer because of character restrictions for comments. Also, I am not a statistican but a molecular biologist who did quite some reading/courses on stats.

Two scRNA-seq tools, namely SingleR and RCA (Reference Component Analysis, Li et. al, 2017) compute correlations of bulk transcriptomic data to single cells; for this they use Spearman correlation and Pearson correlation respectively (The two tools use the same bulk transcriptome reference). I am using these two tools for different purposes so I did not really compare their results, however, single cell clusters that I compute using the reference components (calcuated with RCA) is in good agreement with single cells labels computed with RCA (adjusted rand index of 0.81).

In addition to that, please take a look at scMatch: a single-cell gene expression profile annotation tool using reference datasets, they claim (scMatch computes bulk - single cell correlations like SingleR and RCA):

  From this comparison we found Spearman’s correlation coefficient
  outperformed Pearson’s correlation for all six libraries.


Identify non-coding regions from a genome annotation

Getting the non coding regions of a protein coding transcript, sounds like you are looking for UTR.
UTR has its own feature in the gtf file. So you can do this:
$ awk -v FS="\t" '$3=="UTR"' gencode.gtf

If the gtf file is compressed use this instead:
$ zcat gencode.gtf.gz | awk -v FS="\t" '$3=="UTR"'

BTW: Why are you using such an old release of gencode? The current version is v29.
As @Wouter tells you, the non coding region of a genome is the complement of the coding regions. Coding regions have its own feature in the gtf file. You can get them with $ awk -v FS="\t" '$3=="CDS"' gencode.gtf. Reading the manual for bedtools complement is your task.


I cannot comment on your specific manuscript since I obviously cannot read it but in general the question is whether a study does add an interesting novelty.
The problem with pure in silico studies is that there is essentially no confirmation towards the in vivo relevance. Drug companies simulate and also test many thousands of compounds in a high-throughput fashion on a daily basis so it is often questionable to what extend individual pure in silico studies without in vivo confirmation do actually add to the bigger picture.
That might be the reason your work was rejected but as said this is a rather general comment and I can neither judge content nor quality of your work so please do not feel offended and read it with a grain of salt.


Why is SMALT better for microbial genomics than other mappers?

I have used smalt a couple of times. Smalt is more flexible than bwa/bowtie and can be tuned to be more sensitive to divergent hits, which is useful to certain applications. I heard Sanger did evaluate a few mappers for some non-typical applications (not sure what). They found smalt to perform better. Also smalt was developed at Sanger. I guess they could modify smalt based on their needs, though I don't know the details.


These tutorials on Seurat multimodal data and the wrapper Seurat data are easy ways to start.

The wrapper has some cite-seq data preinstalled making it easy to work with benchmarked data sets

If you are using Seurat u can just as well have a look at Signac to start working with some sc-atac seq


It depends on the tool and the use case but generally they should be.
It's not as important with short reads as only a small portion of reads have supplementary mapping (unless you're looking at breast cancer post chromothripsis, for example). But imagine the case where you're working with long reads (as I do). I have samples with a mean read length of 60kb. The human reference is structurally different to any given sample meaning that the chance of any given long read having supplementary mapping is high. If supplementary mapping didn't count towards allele counts then many would be artificially low. This is my use case, others are different.
I would count them unless you have some reason to suspect that they are are artefact of the editing and or the mapping. Sorry to never give a straight answer but that's bioinformatics! From what you described I would guess 99% are real but there's normally some noise and if that's of concern your code should include some QC. I.e, consider if multiple reads support the sup mapping and their quality.


levels of factors in the design have non-unique level names after make.names() is applied

There's a good reason to not allow + and - in the columns. It's because once you do that you get coefficients that can then start with a + or -, which is illegal in R and most other languages. In R, legal variable names consist of:

  A syntactically valid name consists of letters, numbers and the
       dot or underline characters and starts with a letter or the dot
       not followed by a number.

That doesn't really explain the error message you saw, of course. That can be solved by thinking about what R does to invalid character in names. Namely, it replaces them by a dot. So + and - get converted to X. and X., which are of course identical and are what's really leading to the error message your saw. If you had used +cre and - it would have worked.

As an aside, you'll want to think in terms of positive and negative Cre if you talk to a biologist, since that's the appropriate terminology for it.


Save output to a specific folder and/or with a specific prefix in Cancer Genomics Cloud

If the output is File, you can add this step https://github.com/Barski-lab/workflows/blob/master/tools/rename.cwl to rename the output file.


Snakemake: Migrating from deprecated cluster.json to new profiles.yaml

It looks like you can just bundle your cluster config in its own YAML file, and reference that from within a profile YAML.  Here's a clumsy little example using a fake job submission script: https://github.com/ressy/example-snakemake-profiles

The current documentation makes it sound like cluster configuration support is going away and being replaced by profiles:

  While still being possible, cluster configuration has been deprecated by the introduction of Profiles.

...but from those example profiles, that might be misleading; here's an example template for a profile YAML that references a separate cluster-config YAML:

https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile_name%7D%7D/config.yaml

...
cluster-config: "$((INSTALDIR))/overwrite_cluster_param.yaml" #abs path
...

That cluster-config YAML has the same default and per-rule entries just your JSON version does:

https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile_name%7D%7D/overwrite_cluster_param.yaml

...
# specific parameters for certain rules, which need more time/memory

#run_assembler:
#  queue: bigmem
#  time: 100
#   threads and memory definen in config file

So I think you can just set up a pair of YAML files accordingly including options for your rsem_quantify and star_map_human rules.  From my own test it looks like it handles the cluster argument templating just as it did before when using a profile.  (If this is correct we should nag them to make the new documentation more explicit though.)


RNAseq biological replicates not clustering in PCA plots

PC1 is 81% of the variance?
This PCA plot confirms that different tissues are different.  You probably already knew that.  I'd make more PCA plots that are tissue specific.  That will be more informative than these.
Personally, I'd also not do DESeq on all of these samples together, unless the goal of your experiment really is to learn the differences between these two tissues.


Looking for a good tool for repeat masking in long sequences

RepeatMasker can downloaded an run locally. There's then no size limit, though do note that it's quite slow for mammalian-sized sequences (you'll want to split it into chunks and run in parallel in such cases). For a small 500kb sequence that won't pose an issue.


Time and memory efficient processing of many intervals for a bigWig file

Assuming you can make a BED file of introns, you can then use the pyBigWig module in python:

import pyBigWig
import numpy

bw = pyBigWig.open("some file.bw")
bed = open("introns.bed")
for line in bed:
    cols = line.strip().split("\t")
    vals = bw.values(cols[0], int(cols[1]), int(cols[2]), numpy=True)
    effLen = cols[2] - cols[1] - (vals < 1.0).sum()
    # Do something with this.

One can make that parallel (with the deeptools API) should it prove too slow. This will only use an appreciable amount of memory for very large introns, since there you need to store the value at each base. One could make this much more memory efficient by using bw.intervals() instead, but that'd require quite a bit more code (you'd need to determine the number of bases overlap per interval and then assume that either each base you want to query has an overlapping interval or, if not, keep track of that fact).


Using a simulator for example MinION output doesn't make sense when there are an abundance of freely-available FASTQ and FAST5 sequences for general purpose use (and especially for the purposes of testing out algorithms). Simulation of sequences can only ever express the components of DNA that are explicitly modelled, which is never going to fully represent reality.
As long as you're doing your own private research, you can use the FAST5 sequences from Clive Brown's own genome:
https://github.com/nanoporetech/ONT-HG1
Bear in mind the license attached to that data:
https://github.com/nanoporetech/ONT-HG1/blob/master/LICENSE.md

License to use for research use only (i.e. non-commercial) governed by the law of England and Wales.
You agree to this License by downloading the data on this site whether in whole or in part.
You cannot re-distribute modified copies of these data, or derivative works, without my express written permission. Where such permission is granted, you must cite the original work and its Copyrights and you must include this license with any work that includes any of my genome.

For a more open license, you can use the nanopore-wgs reads here (derived from human cell line NA12878):
https://github.com/nanopore-wgs-consortium/NA12878

We encourage the reuse of this data in your own analysis and publications which is released under the Creative Commons CC-BY license. Therefore we would be grateful if you would cite the reference below if you do.

If you don't want human data, there's mouse parasite data here:
http://www.ebi.ac.uk/ena/data/view/PRJEB20824
And E. coli data here:
https://www.ebi.ac.uk/ena/data/view/PRJEB22772
... and if you just want 8h of sensor noise as a FAST5 file, I've got that here in one of my DropBox folders (as a ~70MB file):
https://www.dropbox.com/s/20j1anlz0motrs4/MINICOL233_20170118_FN__MN15625_sequencing_throughput_ONLL3127_84599_ch320_read10302_strand.fast5?dl=0
Reminder: if you do ever present the results of analyses based on a data source, it's a good idea to acknowledge the source (even if not expected). If nothing else, it will let the creator know that someone cares about their data, and perhaps encourage them to release more in the future.


Chromosomer only exists for Python 2. You should thus be able to install it via

pip2 install chromosomer

But:

Your Python installation is a bit screwed up: there are Python 2 and Python 3, which are unfortunately incompatible.

On your system python and pip seem to be aliases for Python 3, which I’d strongly recommend against (lots of tools will break). Instead, they should alias Python 2. Python 3 should be invoked via python3/pip3.


How to simulate "base error rate" in art_illumina?

You need to specify an error profile to do this. However, you can shift quality scores of an existing profiles to get 10% error, which is easier than defining a complete profile.

For example, Illumina HiSeq 1000 has the error rate of 0.1% and you want 100 time more errors, so you should shift down the quality scores by 20 by setting both qShift and qShift2 to -20. The following code will simulate paired end reads with 10% error.

art_illumina -ss HS10 -i reference.fa -qs -20 -qs2 -20 -sam -p -l 100 -f 20 -m 200 -s 10 -o output


I have 23andme text files and would like to convert to SAM/BAM format

Having done 23andme myself I can tell you that your variant file, which contains SNP genotypes, cannot be converted to a bam file. It does not contain the same information as a bam file. It may be helpful to familiarize yourself with those filetypes and the technologies used to obtain them. A SAM/BAM file contains alignments of reads obtained by sequencing. Your standard* 23andme analysis doesn't have reads, as it uses a SNP array for determining the alleles at predefined positions. Markers which are not included in your 23andme file cannot be added by other means because those were not genotyped by the array. Therefore also Y-DNA cannot use your 23andme file or add additional markers.

*For the sake of completeness: I am aware that under rare circumstances 23andme can perform WGS for additional variant typing, but I do not assume this is the case for OP.


Kallisto psudoaligner quant not giving me expected output

As you surmised, kallisto requires a transcriptome fasta file, not a genome fasta file.


remove sequences from fasta file matching a string in the header

you may use simple R commands, using the 'seqinr' package

#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("filename.fa")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), ";"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c("species_name1", "species_name2", "species_name4")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fa")


How to generate a consensus sequence from a multi-reference BAM file?

This post on biostars has several suggestions.
THe one that seems most straightforward, starting where you are, is to use vcftools:
cat ref.fa | vcf-consensus file.vcf.gz > out.fa

Based on that I think that you should also be able to solve this problem using alfred and GATK.
If all of these suffer the same problem of needing a single sequence (which I very much doubt that they do), you could split the bam / ref files into multiple files, run whatever consensus workflow, and then combine them again. But that doesn't seem like it should be necessary.
See also this post and this post.


multiBamSummary output chromosome region annotate issue

Assuming you used a BED file as input then the output is in the same order, so just merge the two files together.


You don't have any continuous fill scales, since you're filling according to Condition, which is a factor. This is leading to the Error: Discrete value supplied to continuous scale error. You just need to combine your previous scale_fill_manual(values=...) command with the labels= part of your incorrect scale_fill_continuous() command:

ggplot(dataset, aes( x = Norm, fill = Condition, y=log2(measure))) +
    geom_boxplot() +
    scale_fill_manual(values=c("brown1","darkolivegreen4","burlywood3"),
                      labels = paste("condition1", "condition2", "condition3")) +
    theme_classic() +
    labs(title="Comparison",x="Norm", y = "counts")

Something along those lines should work.


What are common ways to calculate gene length for TPM calculations?

I think the most common way is to use something like Kallisto which figures out what proportion of reads are likely assigned to each transcript variant, and calculated TPM  based on that.


It's not clear exactly what your goal is. It might help to have an idea of how you are trying to use this data; and it's probably not a bioinformatics question, as ethnicities are social/political rather than biological constructs.

However, there exist some controlled vocabularies used by LoC for ethnic groups; that looks a bit old so it might not be the best source.

I am by no means an anthropologist but it does look like there are anthropological databases devoted to e.g. ethnography that may have complementary information.

Wikipedia has a simple list of ethnicities that notes, instructively, "There has been constant debate over the classification of ethnic groups."

I think that you are unlikely to find a single authoritative answer that is safe to use as a source of truth. I would normally suggest to try another stackexchange site, though it looks like there is none for sociology or anthropology. There is no anthropology tag available on this site, but you might try posting in e.g. history.stackexchange.com and adding an anthro tag.


Do any publicly available databases detail protein structure and functional domains?

Some of this information (at least some domains, active sites, etc) is available from UniProt.

If you want to download their whole database, you can search without specifying any terms and then click the Download button.


Change dataframe to load all the files into a list:

wd <- "C:/all_your_dataframes_in_one_folder"

setwd(wd)

Then you get the list of files you want from your folder

files <- list.files()

Put into a list with a loop

dflist <- list()
for(i in 1:length(files)){
 dflist[[i]] <- read.csv(files[i])
}

Then use tidyverse solution

library(tidyverse)

map(dflist, ~.$GeneID) %>% #creates a list just of the column of interest
    reduce(intersect)  #applies the intersect function cumulatively to the list


How to filter the PDB for single-domain proteins with full length 3D-structures solved

If you want to reproduce the model in that paper, you should email the corresponding author asking for the exact set of PDB IDs/chain IDs/residue ranges they used. There is too much flexibility in the criteria as given to get the exact set any other way.

If you want to assemble a similar set for your own work, the RCSB PDB advanced search is very powerful and would let you filter e.g. by max sequence length or for only monomers.

However the PDB doesn't deal with domains (regions that evolve and fold independently) directly. A PDB chain may contain part of a domain, one domain, or multiple domains.

If you want a dataset of protein domains you could use CATH, SCOP or ECOD.

As a side issue, please remember to remove overlap of homologous proteins between the training and test sets. People entering the field often neglect to do that. See for example this comment piece.


How do I obtain a vcf read from a pyvcf.Reader?

I think vcf_reader.fetch is what you want.  You have to index the file first though.
If starting from regular text VCF you can run:
bgzip example.vcf
tabix example.vcf.gz

Then the Python code is almost the same as yours:
import vcf
vcf_reader = vcf.Reader(filename="example.vcf.gz")
vcf_read = vcf_reader.fetch("NC_023663.1", 7803980, 7804000)
for record in vcf_read:
    print(record)

For this example VCF you get:
Record(CHROM=NC_023663.1, POS=7803981, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803982, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803983, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803984, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803985, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803986, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803987, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803988, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803989, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803990, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803991, REF=A, ALT=[G, <*>])
Record(CHROM=NC_023663.1, POS=7803992, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803993, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803994, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803995, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803996, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803997, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803998, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803999, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7804000, REF=T, ALT=[<*>])

API reference:
https://pyvcf.readthedocs.io/en/latest/API.html#vcf.Reader.fetch
Note the "zero-based, half-open coordinate system" it describes (which is why 7803980 in the second Python argument becomes POS=7803981 from the VCF but then it stops at POS=7804000).


Ok I think I roughly get what you want to get. Below I simulate something that's like your data because you provided two identical matrices.

#3 TE, 10 tissue, random values
# number of tissue
M1 = matrix(rnorm(30),3,10)
rownames(M1) = paste("TE",1:3,sep="_")
colnames(M1) = paste("Tissue",letters[1:10],sep="_")
## make gene matrix, , random values
M2 = matrix(rnorm(30),3,10)
rownames(M2) = paste("Gene",1:3,sep="_")
colnames(M2) = paste("Tissue",letters[1:10],sep="_")
## here we make TE1 highly correlated to Gene1
M2[1,] = M1[1,] + rnorm(10,0,0.3)
## here we make TE2 highly anti-correlated to Gene2
M2[2,] = -1*M1[2,] + rnorm(10,0,0.3)
#here we calculate for each TE, its correlation with each gene.
#credit to @haci
COR = apply(t(M1),t(M2))
# you can look at the correlation
COR

What you need is the last three rows. Remember, you can only calculate correlation between the gene and the TE, and you do this across all tissue. For example, if gene1 and TE1 has a high correlation, this tells you TE1 expression can predict gene1 expression in most tissue.

It's meaningless to calculate correlation between one data point and another. You need more than 1 value. For example, if you do cor(10,5), it doesn't make sense.

Hope this is what you are trying to do..


Comparing gene abundances between metagenomes

I would avoid using assemblies to answer this question, as there's no guarantee that you will be able to assemble your genes of interest; you can however estimate their abundance even if they are relatively rare.
How I understand your question as being one of estimating the abundance of either some specific genes (e.g. butyrate metabolism genes) or all genes in a microbial community across multiple samples for comparative purposes. In other words, not 16S or marker gene analysis for the purposes of estimating organismal abundance, which is a rather different problem (though in that case I would still not use an assembly).
A more standard workflow is:

align metagenomic reads against some existing database of genes annotated appropriately.
estimate the number of reads aligning against some gene or ortholog in some way (using for example KEGG Orthology or similar).
use counts from (2) as input to some statistical procedure, possibly summarizing across functional categories.

Some examples of how this has been done are here, here, here. I am sure that there are more recent/relevant references but I haven't been following the field closely in the last few years.


Scaling by linear regression against the number of reads

I don't know if this question has been solved already, but what they try to do is equalize the depth of sequencing for each cell. Therefore, they scale for the total number of reads. If you regress out (via linear or negative binomial regression) the differences in the number of reads per cell, you end up with cells that have been sequenced with the same depth.

In my opinion they use the ribosomal genes in the same way. They are considered to be a kind of housekeeping genes, which you can use to equalize the sequencing depth.

And yes, the log in this context is log2, used to get fold-change values instead of counts.


After a bit of digging, finally found exactly what I was looking for. I'll share here in case someone else was looking for the same thing.

I used a standard Nucleotide BLAST (blastn), coupled with MUMmer. Short tutorial.

Assuming that xn & yn are the couple of scaffolds I want to check if they are assembled together in the Pacbio assembly (a multi-fasta file):

1) first, I blasted the pair of Illumina scaffolds (as separate queries) to the PacBio genome (as subject):

blastn -query ../illumina_asm/scaffold_x1.fa -subject ../pacbio_asm/genome.fa > out_scaffold_x1; blastn -query ../illumina_asm/scaffold_y1.fa -subject ../pacbio_asm/genome.fa > out_scaffold_y1

The candidates PacBio scaffolds shall appear in the top blastn hits in both out_scaffold_x1 & out_scaffold_y1.

Let's assume we found one common scaffold only - corresponding to the fasta entry pacbio_scaffold in the PacBio genome.

2) extract the sequence of the target scaffold from the PacBio assembly, using samtools faidx:

samtools faidx ../pacbio_asm/genome.fa; samtools faidx ../pacbio_asm/genome.fa 'pacbio_scaffold' > pacbio_scaffold.fa

3) concatenate the Illumina scaffolds to analyze:

cat ../illumina_asm/scaffold_x1.fa ../illumina_asm/scaffold_y1.fa > illumina_scaffolds.fa

4) compare the Illumina scaffolds to the PacBio assembly using nucmer from MUMmer:

nucmer --prefix=pacbio_scaffold pacbio_scaffold.fa illumina_scaffolds.fa

this will generate a pacbio_scaffold.delta file, used in the next step

5) filter the delta file to keep only the most meaningful alignments file using delta-filter from MUMmer (avoiding confusing messy plots later):

delta-filter -m pacbio_scaffold.delta > pacbio_scaffold.delta.m

6) use mummerplot from MUMmer to achieve a graph:

mummerplot --png pacbio_scaffold.delta.m -R pacbio_scaffold.fa -Q illumina_scaffolds.fa --prefix=pacbio_scaffold -large -layout

On the x-axis, the PacBio scaffold they have in common. On the y-axis, the two Illumina scaffolds (separated by the horizontal grey abline). Colors correspond to the orientation.

Example of Illumina scaffolds linked in the PacBio assembly:


How to fix RepeatMasker.lib.nsq missing error when running RepeatModelor?

What that code snippet is doing:
That snippet you found on github (where is that from exactly?) looks like it tries to run RepeatMasker's configure script, feeding it some options that are needed by RepeatModeler.  The cryptic third line generates and dumps a block of text to a file named tmp:

2
path/to/bin/rmblastn

5

And then passes that in to RepeatMasker's own configure script as though it were entered on the keyboard.  It's the equivalent of hitting enter to accept a default, entering "2" to select the RMBlast option, typing in the path to rmblastn, hitting enter to set RMBlast as the default, and then typing "5" to save and quit.
That all looks fine, but I think this might be a problem with your (and mine, see below) RepeatMasker install rather than RepeatModeler itself.  These tools are new to me but I'm interesting in seeing this work too, so I gave a shot at troubleshooting it here.
The way I was able to get it working (note I assume you're working with conda packages here) was:

pushd $CONDA_PREFIX/share/RepeatMasker
./util/buildRMLibFromEMBL.pl Libraries/RepeatMaskerLib.embl > Libraries/RepeatMasker.lib to "automatically create the old RepeatMasker.lib FASTA style database" as the script's comments put it.  NOTE: you might want to first delete or move the old empty file rather than alter it because of how conda shares data between environments.  If you update it in place (which > will still do), you'll alter the file for all conda environments using the package, and conda may complaint about unexpected file sizes when checking packages.
perl configure to do RepeatMasker's setup including some makeblastdb commands using RepeatMasker.lib.  (Note this is RepeatMasker's configure script and not RepeatModeler's.) You can type in the options manually or you could probably use that code snippet you found.
popd and run RepeatModeler/RepeatClassifier and hopefully it finds the BLAST library files now

All the gritty detail:
I see a conda path in your logs.  Are you using the bioconda packages named repeatmodeler and repeatmasker?  It looks to me like the repeatmasker package in particular includes an empty RepeatMasker.lib file and that's the cause of the trouble.  I don't see any way that RepeatMasker.lib is set up automatically via the configure script, and the BLAST database files like RepeatMasker.lib.nsq that RepeatModeler expects are built from that (empty) file.
The Please rerun the configure program in the RepeatModeler directory message is a red herring, I think, in this case because the $CONDA_PREFIX/share/RepeatMasker/Libraries/RepeatMasker.lib file does exist, just not the other files for a database it assumes should exist.  If I instead run configure for RepeatMasker and update its script to not hide error messages, I see this coming from makeblastdb:
Building a new DB, current time: 10/26/2020 16:10:07
New DB name:   /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib
New DB title:  /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 1000000000B
BLAST options error: File /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib is empty

And it's true, that's an empty file provided by the conda package (the metadata in $CONDA_PREFIX/conda-meta/repeatmasker-4.1.0-pl526_0.json shows size_in_bytes: 0 for that one so that's apparently "correct" for the package.)
I saw a reference to that filename in a buildRMLibFromEMBL.pl script in the RepeatMasker share directory, and running that produces a FASTA file:
./util/buildRMLibFromEMBL.pl Libraries/RepeatMaskerLib.embl > Libraries/RepeatMasker.lib

Then ./configure sees a non-empty .lib file and builds the database, and RepeatClassifier finally sees the resulting database files.
But why is this needed?
Shouldn't RepeatMasker come with a non-empty FASTA file if it's needed by its own configure script?  I can see bioconda's build.sh doesn't do anything in particular to handle that file, but, maybe it should?  This is probably where somebody who knows these tools better can chime in.)


Regressing out unwanted sources of variation in single cell RNA-seq data

Cleaning data before doing analysis can be more important than the method of data analysis. For UMIs, there's an obvious cleaning step that can be done: it would be better to filter out duplicates (prior to generating counts) than to try to incorporate that information in the count analysis.

I'm not familiar with single cell analysis, but I am aware that DESeq2 can be adapted to work with single cells, and incorporating additional variation is just a matter of adding an additional covariate into the DESeq2 formula model.


The direct answer is yes. If you have never worked with hdf5 you can start  here.

But, cellranger outputs several hdf5 files. You can extract the count matrix from the  "/outs/filtered_gene_bc_matrices_h5.h5" file.

They provide a lot of information on how to use the hdf5 format so their their documentation is a good place to start.

In addition: you can also load the data into Seurat, monocle or scanpy and follow their tutorials to analyze the data.


Putting labels of different sizes on one PyMOl Object

'label_size' is a object-state-level setting, which means that you will have to rely on the 'create' command to create a new object for every different label size.

It can certainly be automated, and for this I'd recommend getting into Pyhton Scripting for PyMOL


How to cluster the human genes by pathways/system-biology/metabolic properties?

To measure if two genes are functional similar I developed the BioCor package. It calculates a similarity score between genes by the measuring the amount of pathways  shared. However, it doesn't take into account the strength of the connection between two proteins. You could multiply the similarity by a value between 0 and 1 according to the number of evidences of the relation between the two proteins.

It is optimized for large amounts of data, you could even perform them in parallel, but I doubt you find information about pathways for most of the genes (In my informal test, I think that I could calculate similarities for about 1600 genes not the ~20000 of the human genome).

When I am interested in co-expression data I calculate the co-expression by correlation and add the similarity of the genes calculated with BioCor. The resulting distance matrix can be used for clustering. Maybe the section about using weighted co-expression network analysis and BioCor can be helpful.

To have a "biological pathway clustering of genes for lymphoblastoid or breast cell lines", you would need to have a specific signature/pathway of the cell lines. Currently I am not aware of accepted standard signature database. The single cell atlas could be useful to find expression of your cells, but I don't know if they provide a signature for cell line


Statistical approach to compare the SNP genotyping data among set of individuals

PCA sounds like a good start to see a rough population structure among your 96 individuals, this step is half a sane thing to do, but also it can give you a hint how many distinct populations there are in your dataset. I have not personally run PCA on SNP data, but it seems that there are quite lot of tutorials online. For instance this one looks  nice.

Once you know how many categories you can expect in your population, you can use Structure to classify individuals to the categories.


Short reads vs long reads with regards to break points

A "break point" typically denotes one end of a structural variation, such as an inversion or deletion. You can simply think of it as the position where the normal chromosome was broken to allow the variation. One can find these with short (Illumina) reads relatively easily since with paired end sequencing you'll have pairs spanning these breakpoints with odd relative orientations or insert sizes or such. Note that this is not a particular strength of short reads, since long reads are probably even better at this (they are not as good for finding SNPs at the moment though).


Extracting all reads from bam file which match read IDs in another file

It is still slow but grep has a -f option to take in a file

samtools view inbam.bam | grep -f read_names.txt > read_locs.txt


Relationship between sequencing lane and ngs dataset

If you have used adapters containing barcode sequences (that identify each patient), you would presumably end up with one single tube with your libraries in it right?

Correct, you have to prep each sample into a library.

  If so, you then "put" the contents of that tube along with DNA polymerase and nucleotides into one single lane of a flow cell?

You could, or you could pool it with other samples (assuming they have compatible indices). As @swbarnes2 mentioned, you don't see a lot of samples run on a single lane these days. Sequencing centers like to split things across lanes in case on lane fails for some reason or there are bubbles in the flow cell or other technical issues.

  Your provider would then give you raw FASTQ files that have been identified using this barcode?

Typically, yes, and if your sample was run on multiple lanes they'll often be merged for you automatically. This isn't always the case. I've seen a number of people receive samples either still split by lane (annoying, but you can just cat things together) or, much worse, not demultiplexed at all. In the latter case, people often receive 3 (or 4) fastq files per lane, one for each read and another for each index read. It's rare that sequencing facilities do this, but it happens.

  If you do not use adapters with barcodes, does that mean you would have to use a sequencing lane per patient? If so, that would vastly increase the price you pay your provider, would it not?

You could sample a sequence per lane. Whether that makes economic sense is up to you to decide.


1D² data that is called by Albacore generates two overlapping sets of files. Here's what Devon O'Rourke mentioned on the ONT community forums and on Twitter:

  Following Albacore basecalls on a 1d2 library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D2 script. I'd like to use the reads generated from these scripts in a genome assembly.

  After a little bit of grep searching, it seems like there are overlapping reads among the 1D and 1D2 directories - namely, those 1d reads that generated 1d2 reads!

Link here, for those who have ONT community access.

I've recommended generating a list of the 1D^2 reads, then excluding those reads from the 1D fastq files. I have a couple of scripts that can help for this:

$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

More information on this can also be found in this StackExchange answer.


Single-cell RNA sequencing (scRNA-seq): filtering cells by transcript counts, how to choose cutoffs?

By reading this tutorial:

  Wells with few reads/molecules are likely to have been broken or failed to capture a cell, and should thus be removed.

This answers the first part of the question.

For the second part, i.e. how to decide cutoff, it still seems a bit arbitrary to me, but some authors trim the cells on the left tail of the highest mode, in order to remove other modes (likely due to noise) on the left tail, considering the distribution of cells dependent on the sum of reads (molecules per cell, first plot).


Aggregate sequencing/mapping/etc. metrics from cellranger across Illumina samples

I wrote a simple script to do it. Not sure whether that's what you want.
#!/bin/bash
#CellrangerSummary.sh Sample1 Sample2
csvtool transpose $1/outs/*summary.csv | csvcut -c 1 > summary.csv
for sample in "$@"
    do
        csvtool transpose $sample/outs/*summary.csv > tempfilesforcsvtool.$sample
        csvjoin -c 1 summary.csv tempfilesforcsvtool.$sample > temp2
    mv temp2 summary.csv
rm tempfilesforcsvtool.$sample
    done

csvformat -T summary.csv > summary.tsv
rm summary.csv
```


Extracting expression data in r from GEO and got a S4 object

Well, if you use a SOFT formatted files:

#You can read your file
gds <- getGEO(filename="GSE9838_family.soft.gz")

# and check it
class(gds)
show(gds)

#then you usually would try to use  an exprs() accessor to retrieve the expression matrix:
exprs (gds[[1]])

# but in this case you will get an error revealing that it is a GSE datastructure
# eg. is composite structure having GSM and GPL object lists:
names(GSMList(gds))
names(GPLList(gds))

# now GPLList is probably not interesting in this case
# but you can look at the GSM-objects from GSMList:
first <- GSMList(gds)[[1]]
first
class(first)
Meta(first)

# at this point if you read the results you´ll notice that
#    title
#    [1] "X_intercept"

# and
#   An object of class "GEODataTable"
#   ****** Column Descriptions ******
#   [1] Column      Description
#   <0 rows> (or 0-length row.names)
#   ****** Data Table ******
#   data frame with 0 columns and 0 rows

So here you would probably wonder if the data you want is these objects. If you go back go GEO page for GSE9838, you can find that "Raw data provided as supplementary file. Processed data is available on Series record". So one may want to check series record instead:

# you can query the GSE9838 like this:
data = getGEO("GSE9838")

# then you can check what you´ve got
names(data)
#[1] "GSE9838_series_matrix.txt.gz"
class(data$GSE9838_series_matrix.txt.gz)

# and check what´s inside:
exprs(data$GSE9838_series_matrix.txt.gz)
# but apparently it just shows one header row.

So it wasn´t there either. And if you check the size of the file, you´ll find that it is just 4.8 kb. My educated guess would be that the expression values are not in the dataset itself but probably in some of the supplementary files?

If you want to download them, you can:

filePaths = getGEOSuppFiles("GSE9838")
filePaths
untar("GSE9838/GSE9838_RAW.tar", exdir = "GSE9838/CEL")
list.files("GSE9838/CEL")

...what will allow you to see all the CEL files with the data.


How to generate Beta diversity boxplot from phyloseq object?

You could convert the phyloseq object to a dataframe and just plot it using base R or ggplot2. I'm not exactly sure what the bray curtis df will be called but you can access it using the @ notation then convert with data.frame e.g.
bray_for_plot <- data.frame(phylo_obj@bray)


As @ATpoint explains in his comment, the mapping quality (MQ) applies to read mappings and, hence, is contained in BAM files, but not VCFs. The MQ values carry important information that ideally the variant caller should use when determining the variant quality (QUAL in the VCF) of individual variant loci.
In the LoFreq documentation it says that

Choice of mapper: Bowtie and BWA-SW
BWA-SW assigns very low mapping qualities to mapped reads. LoFreq incorporates mapping quality by default into its variant calling model, which doesn’t help in this scenario and will reduce its sensitivity. Consider disabling the use of mapping quality for lofreq call with -N.

So LoFreq is already aware of the MQ values and takes them into consideration during variant calling.
Since you are using BWA-MEM (and not BWA-SW) you should probably not be concerned about the issue they mention and trust their quality filtering.
If you are sceptic about the quality of mappings, though, you could still filter them out directly after mapping, e.g. using samtools view -q 20 ...


LRT or LRT-like test on cyclical (Sleep) data

I think the problem is that you have time as a linear value rather than a factor. While this naively makes sense (after all, you have to go through 3 hours awake before you can get to 9 hours awake), what you end up doing is fitting a single coefficient to each gene of "linearly changes with time". What you want instead is factor, so you can have a system where things like per and clk can cycle:

time = factor(c(3,9,15,21))

The rest is then correct. Note that LRTs are not asking, "is there a change at any time point?", since there will likely be DE genes according to this that are not DE at any single time point. Rather an LRT such as this is asking, "What are genes that vary significantly over time?" There's a subtle different between those two questions. The former question is actually answered by a bunch of pairwise comparisons (taking the union of the DE genes).


Duplicate genes with RSEM counts: Which one to choose?

There is no one-to-one mapping of gene ids from one database to the other. Ensembl (who maintain Ensembl ENSG IDs), ncbi (who maintain EntrezGene IDs and RefSeq transcript ids) and HUGO (who maintain gene symbols/names) all have different ideas about what a gene is, which part of the genome belongs to which gene etc.

In fact, it's worse than that because HUGO gene IDs refer to the abstract idea of a gene, rather than a particular part of any chromosome.

My suggestion would be to leave the IDs as Ensembl IDs (as these actually refer to a physical part of a chromosome, which is what you are quantifying), right up until you have to report the results to someone who really NEEDS gene names. If you still have duplicates, you could sum, or pick the most interesting, but I would suggest just leaving the duplicates in there, and treat it as thought there are two genes with the same name.


How to filter ribosomal RNA from scRNA-seq data

In the paper mentioned, we used the ScaleData function in Seurat to regress out the number of reads, Rn45s abundance, and percent ribosomal gene transcripts. Ribosomal genes were found with the regular expression ^Rp[sl][[:digit:]].

tiss <- ScaleData(object = tiss, vars.to.regress = c("nReads", "percent.ribo","Rn45s"))

Here's a fuller notebook, and we'll have a better organized repository soon.


Protein ligand docking: how to convert .pdb to .maps.fld?

There are a few steps in between pdb and .maps.fld.
Here is the list of several scripts that can do tasks for you that you downloaded with autodock MGLTools: http://autodock.scripps.edu/faqs-help/faq/where-can-i-find-the-python-scripts-for-preparing-and-analysing-autodock-dockings. Look at the prepare_ files. Also note that the scripts come with a pre-bundled python, but you can install these scripts with conda.
Here are the basic order of things:
prepare_receptor4.py -r protein.pdb
prepare_ligand4.py -l ligand.mol2
prepare_gpf4.py -l ligand.pdbqt -r protein.pdbqt -y
autogrid4 -p protein.gpf
prepare_dpf4.py -l ligand.pdbqt -r protein.pdbqt

partial charge
Your PDB needs partial charges so first you convert it to PDBQT file.
This is true for your protein and your ligand. That is no complicated parameterisation for your ligand —which is not necessarily a good thing!
In this step it is essential to correct protonation if absent each have their own flags.
grid
The you make a grid file (gpf), which contains your grid. To make a box that spans the whole protein simply use the unit cell dimensions form the PDB. For an mmCIF, if you open the header dictionary you need _pdbx_struct_oper_list.matrix using MMCIF2Dict from biopython's PDB submodule. For a PDB file, just search for the CRYST1 line (cf. format).
Once you have have the grid you can run autogrid. This step will create the map.fld
Why a big box is bad for a VS
However, all those who take this approach opt for a box that is huge: as there is no solvent having a huge bounding has little penalty. However, this is considered a very bad strategy for a virtual screen. This not because it uses up computer resources as these are dirty cheap. Structural knowledge is important and targeting the active site of an enzyme stops it, while binding to the surface does nothing unless it's an interface. So manually reading what the stuff are and choosing the box wisely will save you a lot of time downstream in analysis. There is after all the saying "one week in the lab will save you an hour in the library"...
Protein and ligand
Lastly, merge the protein with the small molecule in a dpf file.
Caveat against docking with models
It is not a good idea dock against models. Swissmodel (unless virtually identical), I-Tasser, Phyre or EVFold models. Docking is very sensitive to small structural difference that may have resulted from the modelling, so docking to models is highly discouraged. For Coronavirus protein, dock against SARS, which is highly similar. Or at the very least pay strong attention to the I-Tasser C-scores and discard the bottom 2/3 of their models.
Note about coronavirus
Your objective of docking coronavirus may be a bit too late however.
I did a wee interactive summary of the literature about the solved structure which is the protease. Proteases are very easy to rationally design drugs for.
There are already:

The Covid Moonshot project is an open collaboration between a fragment-screening X-ray facility (Diamond XChem) and many researchers from various disciplines (myself included).
a paper that has an empirically validated coronavirus specific ligand
many many manuscripts about virtual screens in bioarxiv/chemarxiv.
licensed HIV protease inhibitors lopinavir and ritonavir that were initially suspected to be effective —a clinical trial

Also you may want to check out Galaxy project, a tool to run pipelines for genomics, but have started moving into biochemistry from genomics. Eg. https://covid19.galaxyproject.org/cheminformatics/.


How do I efficiently perform a metagenome screen of “all” species?

As you point out, BLAST is not going to work well for reads from a typical RNASeq experiment.

From an RNASeq dataset, you might be able to get away with just a ribosomal screen using Silva, because the ribosome is the most abundantly transcribed gene. This will usually be the case even if polyA selection or ribosomal depletion was carried out on a sample.

However, bear in mind that there will be a lot of false results due to database incompleteness. I used Kraken when I did something like this a few years ago, and ended up with one sequence annotated as a fruit bat (where my subject was a water-living planarian flatworm). Centrifuge would be a better option now; the authors even claim that they can store the entirety of nr in a 70 GB index with Centrifuge (which would be better / more comprehensive than just a ribosomal screen).


Set new Idents based on gene expression in Seurat and mix n match identities to compare using FindAllMarkers

For the first question, you can use ifelse() to create a new column in the meta.data slot:

All_Samples@meta.data$new_column <-
  ifelse(
    rownames(All_Samples@meta.data) %in% colnames(All_Samples_GeneA_Pos),
    "GeneA_Pos",
    "GeneB_Pos"
  )

colnames(seurat_object) provides a vector of cell names in a given Seurat object. Here whatever cell that is in the All_Samples_GeneA_Pos object would be GeneA_Pos and whatever is not GeneB_Pos. To better control the behavior, you can use a "nested" ifelse(); you can have another ifelse() instead of the "GeneB_Pos" bit above.

For the second question, the rationale is the same as above, you can define an extra column with the "identities of interest", then use these identities with SetIdent() and then in turn using two identities of interest in the differential expression call.


Biohackers Netflix - DNA to binary and video

I was quite curious when I saw that as well. I've spent more time than I'd care to admit trying the sort of things you have. I've gotten it decoded now, but I can't really claim any sort of victory, as the problem was pre-solved.
After struggling through some of the same experiments you did, I decided to take a closer look at their explainer video here: https://youtu.be/DMYgjOHgHxc
First I tried to decipher the blackboard sketch around 1:30, but it was quite vague. A closer investigation into the speaker, however, yielded some luck. Googling Dr. Reinhard Heckel brought me to his website here: http://www.reinhardheckel.com/ which shows his most recent publication – a paper on encoding digital data in DNA.
The encoding is relatively complex (something nearly impossible to stumble upon accidentally, I'd reckon) but fascinating. The fragments are indexed and have two, layered error-correcting codes. Perhaps making our job harder (but serving the practical purpose of minimizing homopolymers and unwanted annealing) is the fact that the data is XORed with a pseudo-random noise that shuffles the data. The paper obviously has all of the details.
From my skim, however, I stumbled across this Github repo linked in the paper: https://github.com/reinhardh/dna_rs_coding
As it happens, the README was just updated in this last week, describing how to decode the episode from the file provided. If you have Docker, you just need to copy-paste some commands.
The product is indeed the full and final episode – coming in at a dainty 63.1MB. How did they get 40 minutes of video that compressed? Well, in short, 720x360 at 24 fps. It looks a little awful and is in German only, but is certainly a cool little easter egg. If nothing else, I have a cool paper to read.


Programmatically retrieve graph features from Cytoscape

I was able to extract the graph features that I wanted just using networkx, without going through the intermediate step of Cytoscape.

EDIT: Initially, I planned to import the data into Cytoscape and extract graph features from Cytoscape. However, I noticed that a much easier solution would be to just use the methods of the networkx library. I initially had a few graph features in mind that I wanted to extract like k-core value of each node, and I was able to do this using methods from networkX. To calculate k-core for instance, I used the k_core method. For more information on networkX methods, check the docs


Removing genes with less than a correlation cut-off between two matrices

Nice data and good question.

It is called "outlier" analysis or the analysis of the residual. Its an easy analysis, but . You calculated the standard deviation using the regression slope as the 'mean' and retain data within 1.96 standard deviations. You can plot the data like this and if you are lucky will see the residual forming a nice 'normal distribution'. Outliers outside 1.96 sd can be formally excluded.

Calling (and plotting) the residual from the regression slope will be possible in R, but I use external packages because its easy. BTW you should also try Spearman's coefficient to check you are getting the same answer.


R studio failed to install the latest development version of hdf5r from Github (win10)

Try running this:

source("https://bioconductor.org/biocLite.R")
biocLite("hdf5r")

Just a different method of installing packages besides github. It is more centred around bioinformatics packages as well so it is pretty useful.


What is a sensible number of gene/observations to explain PCA variance?

The general goal of PCA in RNA-seq can be stated as, "I'd like a low-dimension representation of my data to allow easy assessment of the gross structure of my samples, specifically for assessing missing batch effects, samples swaps, etc." Ideally, this low-dimension representation can fit into a plot or two. In other words, we'd like the first 2 (maybe 3) PCs to explain a reasonable amount of variation. A few things to note from this alone:

The number of differentially expressed genes were not mentioned. They have no relevance here.
There was no mention made of the treatment groups clustering together. Seeing this is NOT a goal of PCA.

As you have noticed, the more input rows (genes, assuming you haven't already transposed it) the less variation explained by the first 2 principal components (PCs). On the flip side, the more rows you input the better PCA is describing the structure of your entire dataset.
Question: Do I really care about the gross structure of the entire dataset?
Answer: No, I just want to see problems that I might need to account for.
So you really don't have to use all of the data in the PCA, just the bits that will allow you to see any issues. In general that'll be the few hundred most-variable genes (after some reasonable transformation so you're not just looking at the most highly expressed genes!). The exact number shouldn't really matter very much, something in the range of 300-500 most variable genes should suffice. A common practice is to just randomly pick a nice round number somewhere in or around this range. This is not a number that needs optimization unless perhaps you have a very large number of samples and therefore PC1 and PC2 only account for a small (say 10-20%) percentage of the variation, in which case you're better off looking at either a number of difference principal components or using a different dimension reduction technique, such as UMAP.


It turned out to be easier than expected. You figure out which dbSNP version they used. dbSNP is the classification system for snps based on the rsid format. Pretty much every GWAS publishes their results with the rsids for the variants, as done above as well (MarkerName). So you download one of the dbSNP files from the official site. In this case, I am using version 142 since this was used by the 1000 genomes dataset.
After this, it is a matter of joining the tables. This can be done easily in R (you need about 15GB free memory), we load the dbSNP database into memory. Since we don't need all the columns, we skip some of them to save memory:
#library
library(tidyverse)

#cols
vcf_cols = cols_only(
  `#CHROM` = col_character(),
  `POS` = col_number(),
  `ID` = col_character(),
  `REF` = col_character(),
  `ALT` = col_character()
  )

#this is from the 1kg
#ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/other_mapping_resources/
dbsnp142 = readr::read_tsv(
  "ALL_20141222.dbSNP142_human_GRCh38.snps.vcf"),
  comment = "##",
  col_types = vcf_cols
  )

GWAS_EA = read_tsv(GWAS_EA.to10K.txt,
                   col_types = cols(chrpos37 = col_character(),
                                    chrpos38 = col_character()))

#make chrpos
dbsnp142 %<>% mutate(
  chrpos38 = str_c((CHROM %>% str_replace("chr", "")), ":", POS)
)

#count overlap
(GWAS_EA$MarkerName %in% dbsnp142$ID) %>% table2()

#mutate and merge
GWAS_EA = GWAS_EA %>%
  mutate(
    chrpos37 = str_c(CHR, ":", POS)
  ) %>%
  left_join(dbsnp142 %>% select(ID, chrpos38), by = c("MarkerName" = "ID"))

#update file on disk
GWAS_EA %>% write_tsv("GWAS_EA.to10K.txt")

The two new columns, chrpos37 and chrpos38, have the chr:pos variant names in both versions.


How are PDO and PDX used in computational and predicative models for tumour biology?

Potential pitfall!

I'm not sure about predictive models, but you need to be aware of a potential pitfall in blindly aligning PDX or PDO based sequencing data without first removing contaminating host organism reads, as otherwise these will lead to a lot of false positive variants caused by miss alignment.  In my experience even a small mount of host material can lead to a ten fold increase in called variants due to miss aligned host reads looking like true variants.  I'd recommend using Xenome, source.

Note Xenome was designed for DNA sequencing it's not splice-aware, so for RNA-Seq there is no optimal solution other than removing reads which align to the host organisms genome.  Although the issue here is that in conjunction with a splice-aware read aligner synteny between the host and grafted genomes might create some interesting problems.  However for RNA-Seq provided FACS or similar confirms low-levels of host contamination I expect levels of expression will not be badly affected.  Although this really needs investigating.

Finally and rather annoyingly Xenome produces none-standard FASTQ so you'll need to fix it's output with:

awk '{if (NR % 4 == 1) print \"@\"$0; else if (NR % 4 == 3) print \"+\"$0; else print $0 }' as reported here on seqanswers.


What does final generated models with TM-score of DMPfold signify in "sequence.ema.txt" file?

This is described briefly in the paper. The four quantities are inputs to the estimator of model accuracy:

nres is the sequence length.
neff is an effective sequence count, or how deep the alignment was.
dpscore is how well the final model matches the initial distance prediction.
dpscorenorm is a normalised version of the above.

Don't worry too much about directly interpreting these. They are only given for completion. The important line is Estimated TM-score of ../final_1.pdb is 0.5, indicating that DMPfold just about thinks this is a decent model.

The "sequence.ema.txt" is more raw output giving the probability of the TM-score being in each bin. "rawdistpred.1" is the distance predictions from the first iteration. Always use "rawdistpred.1" as input to the TM-score predictor as this is how it was trained.