If you compare A vs B the genes's fold change will have the opposite sign to B vs A. So will be the gene set up or down-regulated depending on the comparison to take.

The gene set test analyze if a given group of elements is sorted in a certain way in a list (I am talking about the GSE like the one you performed or the one in the Broad Institute). Usually it is used to say that gene set X is up-regulated in the comparison A vs B, that is, gene set X (in general) is more expressed in A than in B. The enrichment is measured in a enrichment score (or a normalized enrichment score, the higher, the clearer trend/distribution group X has. To asses X is really differentially expressed a p-value is also calculated.

Despite that the number of the up-regulated gene sets sum up all the gene sets tested as pointed by juod, the number of gene sets with a a p-value below 5% and 1% is not the same in both comparisons which indicates that you are not really using the same samples in both gene sets.

Check how you performed your comparison. Are the genes in each comparison with the same value and opposite sign?
Check also the software used to run the GSEA analysis (if I recongize the type of report is from a package in Bioconductor [I don't recognize which one now]) as you (or me) might have misunderstood something or there is an error if the comparisons are done properly. For instance it could be that the report of p-values below the threshold would be from only the up-regulated gene sets.


Does the galGal5 chicken assembly have a chromosome 29?

The current ensembl entry doesn't have a 29 either.  The archived ensembl assembly  lacks 29 30, and 31 and 33 and LGE64.

The chromosomes after 30 are tiny, so they might not be visible in a karyotype.  They probably realized that "chr 29" was really attached to some other chromosome.


Here is a sample of code that should get you started.

I created a separate get_abstracts function.  Then that function get called repetitively until less than 100K ids are returned.  If you wanted to get fancy this whole process could be called as a recursive function but that seemed overly complex.

BUFFER_SIZE = 100000

def get_abstracts(start, term):
    """Get N number of records from Entrez
    :returns list: List of ids returned from search
    """
    handle_one = Entrez.esearch(db="pubmed", term=term,
                       retmode="text", retstart=start, retmax=BUFFER_SIZE)
    record = Entrez.read(handle_one)
    idlist = record["IdList"]
    return idlist

list_of_ids = []
start_at = 0

# Iterate the search returning 100000
while true:
    ids = get_abstracts(start_at, 'e coli')

    list_of_ids.extend(ids)

    # Break if less than 100K results were returned
    if len(ids) < BUFFER_SIZE:
        break
    start_at += BUFFER_SIZE

# Do something with list_of_ids


How to calculate module-trait relationship when trait data is in binary format?

Answer from @m, converted from comment:
Personally I'd build a tree. Whats your question patients or the relationship between RNAseq 'signatures'? Its not the definitive analysis, but its good mining.
Pearson is doable but how many Pearson's? You'll end up with a 2x2 matrix. How will you resolve the matrix clustering? A maximum likelihood phylogeny solution would work here making use of the binary format. There are other approaches.
Spearman correlation would be better than Pearson's you end up with a 2x2 matrix.
Answer from @peter-langfelder, converted from comment:
You can use Pearson correlation, it is equivalent to carrying out an equal-variance t-test for the module eigengene between the two groups specified by the binary variable.


I was able to achieve this in the following way:

require(data.table)
cell_names <- vector(mode="character")
for (i in 1:length(neuron_ids)) {
  seurat_subset <- SubsetData(seurat_object, subset.name = neuron_ids[i], accept.low = 0.1)
  genes <- colnames(seurat_subset@data)
  cell_names <- c(cell_names, genes[! genes %chin% cell_names])
}
seurat_subset <- SubsetData(seurat_object, cells.use = cell_names)

Would be interesting to know if Seurat provides such functionality out of the box.


Targeted hybrid capture bam files showing reads outside target regions in my bed file

It strongly depends on exactly what your protocol is. Note that there are a good number of tools for QCing capture data (here, here).
Off-target reads are a standard observation in capture experiments; recall that it is target enrichment, not target perfect purification, as commenter swbarnes2 suggests.
The real question is whether the proportion of those off-target reads is too high for your application, or otherwise suggests problems. Those depend on details of the technique that you will likely have to figure out yourself.
Here are some tables from the first of those QC papers (supplement), showing on-target proportions for various capture samples (sorry, it would have been a bit of effort to get text of the tables rather than images):

These data suggest that on-target reads are 60%-80% of mapped reads. Note that they also lose a lot of reads in mapping, so we cannot exclude the possibility of mis-mapped garbage reads contributing to off-target mapping.


How to merge two files on the first 3 columns?

You can use the same approach as for two fields, just add the third:

awk 'NR==FNR{a[$1,$2,$3]=$4;next} ($1,$2,$3) in a{print $0, a[$1,$2,$3]}' ret lcls


First of all, you aren't actually setting the $outfile variable. Your loop has this basic format:
for var in *; do newVar="something" command $newVar; done

For example:
$ for var in 1 2 3 4; do newVar="new:$var" echo "$var: $newVar"; done
1:
2:
3:
4:

As you can see, newVar is never actually being set. That's because setting a variable is its own operation. So you need to do that, end it with a control character like ;, or &, or a newline, and then run your cutadapt:
$ for var in 1 2 3 4; do newVar="new:$var"; echo "$var: $newVar"; done
1: new:1
2: new:2
3: new:3
4: new:4

Or, with your data:
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

But that would give a weird file name, like foo_S1_L001_read1.fastq_trim.fastq. You probably want names like foo_S1_L001_read1_trim.fastq instead, so you would want:
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

Now, the specific error you show makes me think that isn't the command you ran. Could you perhaps have run this one instead by mistake (note the / instead of \ before the _trim.fastq)?
for infile in /data/HTS_seq/*_read1.fastq; do
 outfile=$infile/_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
-e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

Your error message is:
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

Which means that cutadapt tried to open /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/ as a directory and then look for a file named _trim.fastq inside it. The only way I can think that this would happen is if you ran the command I gave above and not the one you show in your question.


Seurat Merged objects tSNE - How to paint on original IDs?

Below a few lines of code that accompany BC Wang's answer.
After using MergeSeurat the sample name will be added to meta data under orig.ident. this can then be used to color the tSNE either using group.by or pt.shape. The former will show colors for each sample, the latter will color each cluster and give sample id another shape.

path1 <- file.path("path to 10x sample 1")
path2 <- file.path("path to 10x sample 2")

S1 <- Read10X(data.dir = path1)
S2 <- Read10X(data.dir = path2)
S1 <- CreateSeuratObject(raw.data = S1, min.cells = 3, project = "S1")
S2 <- CreateSeuratObject(raw.data = S2, min.cells = 3, project = "S2")

s.combined <- MergeSeurat(object1 = S1, object2 = S2, add.cell.id1 = "S1",
            add.cell.id2 = "S2", project = "2samples")

All data analysis goes here:
RunTSNE etc...

Plot the data

TSNEPlot(s.combined, group.by ="orig.ident")
TSNEPlot(s.combined, pt.shape="orig.ident")


Interpretation of hierarchical clustering with bootstrapping

Good question. You have performed the first step via unsupervised learning and AU bootstraps give great results, orthodox bootstraps give no significant (or very few) clusters. The contrast between these approaches is unusual because they should be somewhere close to each other. If you can resolve that, your analysis is fine. So at present IMO the unsupervised learning needs further investigation.
... deleted

... antibodies ... oh I see. Trees, like you've done, don't work here because of the cross-reactivity creates network phylogenies. This is the antonym of bifurcating trees. If someone has a robust bifurcating phylogeny with antibody data and gets good bootstraps .. utterly amazing. However, if you subject your data to a network analysis you'll find connections all over the place - in practical terms this really messes up bootstraps. Thats the likely reason for you results. There is no point doing a complex bifurcating tree analysis because the antibody cross-reactivity will mess it up. k-means has been implemented as unsupervised learning in these instances and would provide the starting point in using clustering to define a supervised learning problem.


How to get a list of genes corresponding to the list of SNPs (rs ids)?

I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )

A simpler solution, I don't know... but this is at least one Python solution using Biopython's ELink method via NCBI's Entrez E-utils.

The Biopython library is flexible enough (they have an in-depth tutorial worth reading) to modify the code below to to fit project-specific needs.  For example, the mapping between snp id and gene id is lost with this bulk request method.

from Bio import Entrez
Entrez.email = "youremail@gmail.com"

snp_ids = ["rs8179065",
           "rs12899845",
           "rs8179066",
           "rs80357332",
           "rs80356920"]

record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=",".join(snp_ids).replace('rs', ''),
                                  db="gene"))

This returns a record object

[{'LinkSetDb': [{'DbTo': 'gene', 'LinkName': 'snp_gene', 'Link': [{'Id': '374654'}, {'Id': '90381'}, {'Id': '672'}]}], 'ERROR': [], 'LinkSetDbHistory': [], 'DbFrom': 'snp', 'IdList': ['8179065', '8179066', '12899845', '80356920', '80357332']}]

which contains several iterable fields -- 'Link' has a list of the Gene UID's that can be converted to names using ESummary.

for gene_id in record[0]['LinkSetDb'][0]['Link']:
    handle = Entrez.esummary(db="gene", id=gene_id['Id'])
    uid_record = Entrez.read(handle)
    handle.close()

    uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
    print (gene_id['Id'], uid_summary['Name'])
    # ...or something more interesting than print

374654 KIF7
90381 TICRR
672 BRCA1

An option to preserve discrete mapping, is to loop over each snp id and extract the gene(s):

although I don't endorse the method below, especially for long lists of snps, you may flood the server with too many requests, maybe someone more familiar can comment on request limits or best practices for Entrez

for snp_id in snp_ids:
    record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=snp_id.replace('rs',''),
                                  db="gene"))
    results = record[0]['LinkSetDb'][0]['Link']
    for result in results:
        uid = result['Id']
        handle = Entrez.esummary(db="gene", id=uid)
        uid_record = Entrez.read(handle)
        handle.close()
        uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
        gene_name = uid_summary['Name']

        print(snp_id, uid, gene_name)

rs8179065 374654 KIF7
rs12899845 374654 KIF7
rs12899845 90381 TICRR
rs8179066 374654 KIF7
rs80357332 672 BRCA1
rs80356920 672 BRCA1


variant calling on ChIP-seq style data: samtools mpileup with minimal filters

I used this in the past for ChIP-seq data and it generated SNVs:

samtools mpileup \
--uncompressed --max-depth 10000 --min-MQ 20 --ignore-RG --skip-indels \
--fasta-ref ref.fa file.bam \
| bcftools call --consensus-caller \
> out.vcf

This was samtools 1.3 in case that makes a difference.


Why does RATT create a directory with chromosome sequences and how can I stop it?

RATT seems to be a very simple wrapper script that unfortunately does not clean up its temporary files properly. Short of modifying the code yourself, your only option is to perform the cleanup manually after running RATT.

The temporary copies are made because (a) the query is not an exact copy: RATT performs some cleaning. And (b) it needs to extract a FASTA file from the EMBL input for use in downstream tools.

The “mysterious” number behind the file names is simply the process ID of the RATT process. This is appended to make these intermediate files more or less unique (presumably to prevent overwriting existing files).


Multi Omics comparison RNA seq : Condition specific gene signature

From their methods section:

On the basis of these lists of marker genes and receptor–ligand pairs, we inferred potential cell–cell interactions for all pairs of one structural cell type and one haematopoietic immune cell type, quantifying the enrichment for known receptor–ligand pairs among all pairs of marker genes between the structural cell type and the haematopoietic cell type (Fig. 2a). First, we counted all pairs of marker genes for each pair of structural and immune cell types. Second, we calculated the fraction of these gene pairs that were annotated as receptor–ligand pairs. Third, we tested whether this fraction was greater than the fraction of annotated receptor–ligand pairs across all pairs of genes. Fisher’s exact test was used to obtain P values and odds ratios, as implemented by the function ‘fisher.test’ in R. P values were adjusted for multiple testing using the Benjamini–Hochberg method, as implemented the function ‘p.adjust’ in R. Finally, significantly enriched pairs (adjusted P < 0.05) of structural and immune cell types were connected by edges to generate a graph of cell–cell interactions.

In addition to testing for significant enrichment (as described in the previous paragraph), we calculated aggregated receptor and ligand gene signatures (Fig. 2c) based on a manually curated list of immune-related receptors and ligands (Supplementary Table 4). These gene signatures were derived by normalizing expression values for each gene and averaging across all genes in a given set of biologically related receptors and ligands.

The actual values plotted then seem to be from EnrichR, since they mention that elsewhere in the methods.


If you are working under a Linux environment, then you may easily use a piped grep search. For example:

grep deleterious file.vcf |grep damaging >filtered_file.vcf


how to do analysis of a table content SNP?

For working with a relatively small dataset in tabular format, I would recommend dplyr. You can read in a excel file using readxl

library(readxl)
library(dplyr)

df <- read_excel(FILENAME)

To get the number of mutations per gene:

df %>% group_by(`Gene name`) %>%
  summarise(num_mutations=n())

If you want to know if there are mutations in a particular gene:

df %>% filter(`Gene name` == QUERY)


Error in FindIntegrationAnchors- Seurat package

There's a few problems with your code, first, when you do Read10X() it returns you a sparse matrix, and you need to put this into a Seurat object with meta data, before doing the integration.
So for example, i create some example data that is similar to your output from Read10X() :
library(Matrix)
epithelial = Matrix(matrix(rnbinom(100000,mu=200,size=1),ncol=50),sparse=TRUE)
rownames(epithelial) = paste0("gene",1:2000)
colnames(epithelial) = paste0("e",1:50)

habermann = epithelial + rnbinom(100000,mu=20,size=1)
colnames(habermann) = paste0("h",1:50)

Then convert to a seurat object:
data.list=c(epithelial = CreateSeuratObject(epithelial),
            habermann = CreateSeuratObject(habermann))

Then find variable features:
data.list <- lapply(data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Should look like this:
$epithelial
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

$habermann
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

Now perform integration, below I have to reduce k.filter because I have very little cells in this example. In your case, you can simply use the default settings.
features <- SelectIntegrationFeatures(object.list = data.list)
anchors= FindIntegrationAnchors( data.list,max.features = 200,
k.filter=50,k.anchor = 3,verbose = TRUE)


Motivation behind the neighbor-joining distance matrix recomputation

NJ as with all clustering methods there isn't an explanation, It is a hierarchical clustering method but there are many algorithms to cluster a distance matrix. NJ is closely associated with the concept behind UPGMA (unweighted pair group mean algorithm)

The key difference with UPGMA is that after joining two taxa by NJ the reconstructed genetic distance between the node and one taxa versus the other taxa is not assumed to be equal, i.e. the rate of mutation is heterogeneous. All other clustering algorithms I know assume a "global molecular clock" if that makes any sense, the mutation is constant against time. NJ removes that assumption and thats why we like it.


Chromosomal order of supplementary alignment in BAM file "SA" tag

There is nothing in the SAMtags specification that dictates the order of the entries in an SA tag.  It may be entirely dependent on the aligner that you used.

Looking at my own data that was aligned with bwa, I cannot see a definite pattern.  Longer matches tend to be first, but not always.

Your aligner documentation may have more information on how it uses SA tags.


Behavior of `--reference` flag with samtools sort

You would only need to supply a reference when using CRAM output encoding. The CRAM format allows storing entries without explicitly encoding the base calls, since one can instead encode the differences to the reference genome.

With BAM output the --reference option has no effect.


Your dataset was produced on a microarray. The probe IDs in the FASTA file are similar to read IDs produced from other NGS experiments like RNA-Seq. According to the workflow provided by Ensembl the first step is to align your reads to a reference. You'll only keep the reads that match the reference completely. Once the probes are aligned you'll be able to quantify the number of reads mapping to the gene that the probe was designed to attract.

I did find a file providing the database entries associated with each probe ID (one of which is the ID you provided) for the Affymetrix Cotton Genome Array.


Your question can be rewritten as, "how do I get the union of reads in two files?" Since you have already subset the mates, the steps are:

Process each file, storing the list of reads names in each.
Take the union of these.
Iterate through the original/non-subset files and write only the reads with names in the union to the output.

The files will be both subset correctly and synced properly such that bowtie2 or any other aligner can correctly process them.


Is there a way to use BioLemmatizer from Python?

This area is huge in Python because of machine learning text processing, in fact O'Reilly have an entire book on it . However, we don't really deal with words in Bioinfo SE,

You might want to look at the nltk (natural language toolkit) for Python, but for direct parallels with the Java package bio package (your question) I don't think it's something we would know.

You could solve this via subprocessing into Java, and passing the result back to Python. I get the idea that the core ML/ DL used in natural language is set up for Python. Google for example have fantastic tensorflow stack for this.


How do I filter a GFF file by gene type?

grep -E 'gene_type "(processed_pseudogene|lincRNA|unprocessed_pseudogene|antisense_RNA)";' input_file.txt

Explanation

grep -E: use extended regular expression so you don't need to escape things.
This kind of expression (foo|bar|baz) allows grep to include either of these strings. Hence, the total expression will match gene_type "processed_pseudogene";', gene_type "lincRNA";', etc.


pymol script to only select 5 chains within a distance from a reference and save the selection

Firstly, I would preface that the PyMOL command line in the GUI is great, but you can run PyMOL as a regular Python library (calling the appropriate python). This makes life so much simpler for more complex manipulations. You can do really complex manipulations with it. I cannot suggest it enough.
The API version of the commands works only marginally differently, but you can do API calls in the GUI cmd.) and PyMOLWiki documentation is only slightly weaker for them (help() fills the gaps).

But if you really want to use the GUI route...

Let's get the names of the models in the selection

onames = cmd.get_object_list('chains')
print onames

Now, you want to save subsets of 5 chains (technically, objects or models). That is, you do not care about the order they are in and you don't want the same one twice (actually selections with 4 or less), so you want all combinations. Make sure to do the following
    print len(onames)

And google that number as <n> choose 5. That is how many files you will get. I am pretty sure most of these sets will be meaningful and will be useless.
In a script a loop is easy, while in the gui you have to keep it simple or use a list coprehension hackishly.

 from itertools import combinations
 for s in combinations(onames, 5): cmd.save('-'.join(s)+'.pdb', ' or '.join(s))

Obviously, you could alternatively use an enumerate method around the iterator and get an index for the file. But do note f-strings will not work.

 for i, s in enumerate(combinations(onames, 5)): cmd.save(str(i)+'.pdb', ' or '.join(s))


Passing custom distance matrix to T_Coffee


The Retrieve/ID mapping tool from Uniprot allows you to convert from Uniprot to entrez ID and vice-versa. It is available here. Just select the appropriate values in the 2.Selection options panel.


Why does DESeq2 convert numeric columns to factor during differential expression analysis?

Edited:

It's the "as.matrix" that is causing the problem.  When I use that on my sample annotation data, and make dds with that, it turns my numbers into factors


Selection of differential expressed genes

You can usually get away with FDR < 0.1, but that's as high as you can go. This all presumes you're doing follow-up experiments of some sort, of course. I guess you could increase the FDR more, but you're then really increasing the odds that your follow-up experiments will fail. Obviously increasing the FDR will decrease the confidence in the results, it's called "false discovery rate" for a reason.
Yes, fold-change cut-offs are mostly useful when you need to prioritize hits further. Small changes are almost impossible to validate with less-sensitive methods (qPCR) and less likely to be biologically relevant.
You're just going to waste your time if you ignore FDR.

Remember that all thresholds are intended to increase the odds of success in follow-up experiments. If you have infinite time and resources then you can set these to whatever you want, but you have neither so you better prioritize reasonably.

It's likely that you have a large batch effect in your data, effect sizes in cancer tend to be much larger than this.


predict the IC50 after using Schrodinger suite Glide tool for docking

I don't think it is possible just from the docking scores. Docking scores are designed to give you an estimate of the binding affinity (entalphy), but they do not take entropic contributions into account very well. The best you can do are molecular dynamics simulations to calculate the free energy of binding. Which can be computationally quite costly e.g. with GROMACS or CHARMM.
You can estimate the IC50 with molecular dynamics simulations as well, but that is not straightforward. E.g. with Umbrella simulations. An example for a potassium channel and a very simple ligand can be found here.
Another approach would be to build a predictive model. If you have a library of compounds at hand you could use the chemical structures and e.g. the docking scores to estimate the IC50 based on statistics or machine learning. But biology is very messy in general, though if you get a representative dataset with very similar compounds to the one that you have in your dataset then you might be successful. You have to be sure that the binding site and pose are correct though.


Error: Insufficient values in manual scale. 24 needed but only 1 provided?

You must supply genes (or features) as a vector, not individually:

VlnPlot(tumors, c("Cd79a", "Foxp3", "CD14"))

By the way it is probably Cd14, not CD14.


Here's one way:

$ sed -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

The idea is to match a >, then as many non: as possible ([^:]*), and then capture the string between the first :: and the second ::. Finally, you replace the entire line with the captured string.

Once you're satisfied this does what you need, you can edit the original file with -i:

sed -i -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa

Alternatively, set awk's field separator to :: and print the second field only on lines that start with >:

$ awk -F'::' '/^>/ ? $0=">"$2 : 1' file.fa
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*


As the description of the package says:

  The AnnotationHub web resource provides a central location where
  genomic files (e.g., VCF, bed, wig) and other resources from standard
  locations (e.g., UCSC, Ensembl) can be discovered.

If the new resources aren't in the standard locations it won't be listed by AnnotationHub. The easy/common way to have the data accessible through AnnotationHub is make it available through these standard locations (before publications if you want to give early access).

Another method would be to follow this steps of the vignette to add new resources, but it might involve creating a new package.


They are both lentiviruses and share a distant common ancestor.
HIV-1 and HIV-2 are descendents from simian immunodeficiency virus (SIV). The following tree shows the relationships very clearly, from Wertheim and Worobey (2009) Dating the Age of the SIV Lineages That Gave Rise to HIV-1 and HIV-2 PLoS Comp Biol. here. The evolution of HIV-1 is heavily mixed with chimanpzee (cpz) SIV, and there are regions of highly anomalous homology juxtaposed with regions of notably lower homology.


Here's my current understanding of the freemix metric output by verifyBamID. The freemix score is a measure of contamination that is calculated without knowing the genotypes of the individuals. In the paper it is equation (2), and the authors refer to it as the "sequence-only" method. However, in addition to the sequencing reads, it requires some knowledge of SNPs and their allele frequency from the population. They demonstrate in the paper that using the SNPs and allele frequency estimates from a different population actually reduces the estimated level of contamination.

As for intuition for what equation (2) is measuring, imagine what would happen to the number of heterozygous sites detected if DNA from two different individuals was sequenced together. The total number of heterozygous sites would increase because it would count all sites where the two individuals are homozygous for opposite alleles or at least one of the two individuals is heterozygous. This would lead to an increased fraction of heterozygous sites over that expected from Hardy-Weinberg Equilibrium (i.e. 2pq).

Here's the explanation of freemix directly from the author, Hyun Min Kang:

  The key idea of FREEMIX estimate is to use excessive heterozygosity to estimate the level of contamination. Especially for common SNPs, you will observe higher fraction of heterozygous alleles than 2*p*(1-p), and it turns out that you can quantify the contamination very well if you know the population allele frequency already. If you do not have accurate population allele frequency information, than it would be harder to estimate FREEMIX parameters using verifyBamID.

Additional sources:

Verifying Sample Identities - Implementation
Interpreting output files
verifyBamID Google Group


How to programmatically classify a protein according to its genbank feature

Using the UniProt website API (https://www.uniprot.org/help/api) this could be done with a query like
https://www.uniprot.org/uniprot/?query=name%3A%22type%20IV%20secretion%20protein%20Rhs%22&columns=id%2Centry%20name%2Creviewed%2Cprotein%20names%2Cgenes%2Cgo(biological%20process)%2Cgo(molecular%20function)%2Cgo(cellular%20component)&sort=score
or, in tab-separated format:
https://www.uniprot.org/uniprot/?query=name:%22type%20IV%20secretion%20protein%20Rhs%22&format=tab&limit=10&columns=id,entry%20name,reviewed,protein%20names,genes,go(biological%20process),go(molecular%20function),go(cellular%20component)&sort=score
(Remove the "&limit=10" part from the URL to obtain the full set instead of the first 10 entries.)
You can add colums as described here:
https://www.uniprot.org/help/customize,
https://www.youtube.com/watch?v=ado1r8IDm3U (video, at about 1min25)
or
https://insideuniprot.blogspot.com/2015/03
Note that while UniProt tries to improve protein name annotation, the unreviewed section of UniProtKB (i.e. UniProtKB/TrEMBL) is often far from perfect and its accuracy relies on the information provided by the submitter of the nucleotide entry, unless an automatic annotation rule is available to assign a more suitable name (see https://www.uniprot.org/help/protein_names).
Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.


How to import data from cell ranger to R (Seurat)?

Have a look here and the entire tutorial:
https://github.com/hbctraining/scRNA-seq/blob/master/lessons/03_SC_quality_control-setup.md


How many reads has my sequencing run produced on minion?

That's indeed the number of reads and that's quite low. How was your pore occupancy (number of pores sequencing) and your flow cell QC (number of pores good enough for sequencing)? How was your library concentration?

Events is an approximation of nucleotides: it's a guess based on the current signal. As far as I know if you multiply this by 1.8 you get approximately your throughput in nucleotides.

The 4,120,722 is the number of events in the currently selected (mouse hovering?) bin of the histogram.


Yes, Biopython deliberately does not allow duplicate keys with .to_dict(), e.g.
Input file:
$ cat test.fa
>1
ATG
>2
AAA
>1
CG

Code:
from Bio import SeqIO

print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta")))

Traceback:
  File "script.py", line 3, in <module>
    print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta")))
  File "/path/lib/python3.7/site-packages/Bio/SeqIO/__init__.py", line 789, in to_dict
    raise ValueError("Duplicate key '%s'" % key)
ValueError: Duplicate key '1'

However, you can easily write your own function, e.g.:
from Bio import SeqIO

def to_dict_remove_dups(sequences):
    return {record.id: record for record in sequences}

print(to_dict_remove_dups(SeqIO.parse("test.fa", "fasta")))

Output:
{'1': SeqRecord(seq=Seq('CG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[]), '2': SeqRecord(seq=Seq('AAA', SingleLetterAlphabet()), id='2', name='2', description='2', dbxrefs=[])}

Alternatively, you Biopython's .to_dict() accepts a key_function argument, e.g. if you wish to make all the headers unique
from Bio import SeqIO
from itertools import count

c = count()

print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta"),
                    key_function = lambda x: f"{x.id}_seqNum{next(c)}"))

Output:
{'1_seqNum0': SeqRecord(seq=Seq('ATG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[]), '2_seqNum1': SeqRecord(seq=Seq('AAA', SingleLetterAlphabet()), id='2', name='2', description='2', dbxrefs=[]), '1_seqNum2': SeqRecord(seq=Seq('CG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[])}


Smallest group size for differential expression in limma (bulk RNA-Seq)

The problem with most of the methods is that use the gene's variance for each group, which can't be calculated (reliable) when the sample is <= 2. Also, statistically, it would have extremely low power, so the conclusions couldn't be trusted much.

You could also calculate the "raw" fold change by yourself (ie without the variance estimation and adjusting). If you had one more sample per group you could use DESeq2 which was thought for such cases. But from the paper:

  However, if there are two or fewer replicates for a condition, these samples do not contribute to outlier detection, as there are insufficient replicates to determine outlier status.

Which makes me unsure if it will work well for just two samples as it won't help to detect outliers genes.


How to run MaxQuant in command line mode?

Give

mono MaxQuantCmd.exe [OPTIONS]

a try.

University of Oslo description of how to use the tool in their infrastructure.


How to convert the .vcf (imputed) file with GT:GP format to GT:DS?

You can do this in Hail.

Here's the rough code to do it (0.1 versions).

Setup:

from hail import *
hc = HailContext()

Import the .gen file. VCF works too:

dataset = hc.import_gen(
    'src/test/resources/example.gen',
    'src/test/resources/example.sample')

Remap the genotype schema and export to VCF:

dataset.annotate_genotypes_expr('g = {GT: g.call(), DS: g.dosage()}')\
    .export_vcf('/tmp/out.vcf.bgz')

Take a look at the getting started page if you want to try it out!

I should note that you may be able to do QTL analyses in Hail, depending on the method you want to use. See blog post here.


Can open-source software be peer-reviewed and published?

There are whole journals based primarily around publishing open source tools. The primary example of that is "Bioinformatics", where a lot of the open-source tools are published. We've also had luck publishing in the Nucleic Acids Research yearly special webservers issue, since we make Galaxy wrappers around our tools.

You can also publish on tools in venues like Nature Methods, but it's vastly more difficult to get in there if your paper is purely on a tool and it's not something game-changing like salmon or kallisto.

Regarding what's actually needed to publish a tool, you will typically need the following:

Decent documentation. I ask for this as a reviewer and am asked about it when I submit things.
Easy installation. No one wants to spend two hours trying to compile your tool, so make sure there's a conda package or at least a docker container for it.
Test data. Without this your reviewers and users will likely not be happy since they won't be able to try it out on something small.

For the actual submission you'll tend to need comparisons to other tools. I'm generally not a big fan of this, but you will find that some reviewers will demand it (assuming there are other tools that do something similarish to yours).


Using -SP is equivalent to running bwa mem on each of the two mates as if they are single-end reads, but it formats the output as a proper paired-end output, i.e. with all pair-related flags added properly. Without -SP, by default bwa mem forces an alignment of a poorly aligned read if its mate is aligned somewhere. -SP turns off the forced alignment. We use -SP for mapping Hi-C reads.


Why does samtools mpileup sometimes include ref bases (other than ',' or '.')?

Why is 't': 2' in there?

How else do you think the pileup should show you that two reads (in the reverse direction) have a T in that position?


How can I interpret gene expression data from Bioconductor packages?

Processed microarray data is commonly represented as log_2(expression). This data transformation is used because the data more closely fit a normal distribution in log space. With such a transformation, negative values represent very small amounts, rather than the absence of something.

Larger values mean that more stuff is being produced (presumably transcripts), but its significance can only be assessed when combined with other data (e.g. through differential expression analysis).


I cannot access the pipeline from nanopore, but looking at the python script you provided:

with pysam.AlignmentFile(bam_file, "r") as bam:
    for tid in range(0, bam.nreferences):
        ref_name = bam.get_reference_name(tid)
        if len(list(filter(lambda x: x in ref_name, filter_names))) > 0:
            continue

        print(ref_name,
              "0",
              str(bam.get_reference_length(ref_name)),
              sep='\t')

So the script gets the length of the chromosome and by default, outputs the start as 0.Zero based coordinate system is sometimes used in UCSC systems, so I am only guessing the authors of the pipeline use this convention somewhere further down.

As for the coverage you mentioned, that is usually output as a score in the 5th column (bed format). In those situations, usually the coordinates of the alignments in the bam file are converted into bed format. In this pipeline, i think the authors are simply extracting the chromosome information out and putting it into a tabular form.


PDB file downloading: pymol automation vs. manual

To expand on the comment by marcin:

fetch downloads files in mmCIF format by default (https://pymolwiki.org/index.php/Fetch). Not all PDB entries have PDB format files, e.g. due to too many chains. Presumably this is why the change was made, though mmCIF files tend to be larger and hence download slower.
When calling save with the .pdb file extension the structure is converted to PDB format. Perhaps the converter uses the SEQRES records to define the sequence, and hence gives bad files when there are some residues not resolved.

To correct it you could try:

Saving as .cif and reading in with the Biopython mmCIF parser.
Providing the type=pdb argument to fetch, which will force PDB format downloads. Some PDB files may not be available.
Don't do it in PyMol and instead use the Biopython structure downloader (http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc187).

Also, if download speed is limiting you then consider using the MMTF file format, which is a binary format and hence smaller. Biopython can read in and write out MMTF files.


Is the quality score of fastq used somewhere besides trimming/fastqc?

I think almost all alignment programs would use base quality as a factor to calculate mapping quality (See here). These would also be used in Variant Calling to pick reads where high quality bases confidently align as either the reference base or one of the alternate bases.


Besides generic nucleotide aligners, there are also more specialized tools for the alignment of 16S amplicon sequences, e.g. SINA (article, software) which is part of SILVA


Finding simple sequence from reads with significant overlap

Miniasm will probably work for this; it's designed to be an ultra-fast, greedy assembler, which has been used for people doing real-time assemblies from MinION data. Because of its greedy nature, it can misassemble regions where the overlapping fragments are ambiguous, so it's worth using some other method in addition to miniasm to confirm that the assembly is correct.
By default it looks for overlap regions that have a coverage of 3 reads when joining (option -c), but I assume that could be reduced down to 2 reads, as in your example (which has a minimum overlapping coverage of 2).


Bonferroni count = 1 in GO-Slim overrepresentation test (PANTHER) is it possible?

GO-Slim BP, Binomial, Bonferroni

This seems easy you because I suspect have run too many tests at 5%.
If you run 10 tests at 5% rejection threshold, the Bonferroni correction threshold is 0.5% - you divide the number of tests by the rejection threshold.
One possibility is the rejection threshold is given in number of genems: in one instance you've run a pairwise test with several thousand genes so the rejection threshold has been set accordingly. The other test you ran I suspect was a multi-variate test, thus the shift in rejection threshold measured in number of genes.

The altnative is that Bonf = 1 is in fact P = 1, in other words there is no significant difference at all.  If P = 1 then this could be no genes show deviation hence 1770, because 1770/1770 =1
Looking at it again the above seems a better explanation.


Problem with SPAdes installation - no bin folder

Simply unpacking the source code won't actually install much, you'll need to follow its installation instructions. Note that this might be possible on Windows if you use new bash command line, but in general don't expect any bioinformatics tools to work on windows.


Expression of a gene in different groups

You should never use RPKM. It’s simply obsolete in the age of paired-end sequencing, and has been replaced by FPKM (which is, strictly speaking, a synonym).

The linked blog post explains more generally the problems that measures such as FPKM and CPM suffer from. A more robust measure is the TPM (transcripts per million), which scales CPM by the lengths of the transcripts.

Lastly, none of these measures perform very well for comparing gene expression across replicates. Since you have read count data, you should run DESeq2 or edgeR to perform differential expression analysis across all genes, and go from there.

As the comment noted, the paper you cite does not represent scientific best practice. In fact, a competent reviewer would have rejected this analysis.


Why is there a chloride ion in this 3D model?

It is a crystallisation artefact. Namely, the protein are placed in a condition where they fall out of solution without aggregating a gloop (as happens in most well in xstal trial) and to pack as a crystal they need to be placed nice and orderly. It is coordinated by water so is not relevant. A lot of structures have these —DMSO and ions are the most commonly.

Two footnotes.

1ZV7 is a crystal structure, i.e. solved by X-ray diffraction. A model is generally used to described an in silico model, such as those from Swissmodel, phyre, I-Tasser etc. These are generated ab initio (from scratch) or threaded (mapped on top of) against crystal/EM/NMR structures. For Covid I-Tasser has a nice picture of its models.
The best source of info about structures is the PDB. RCSB PDB is more commonly used, while PDBe is more tricker to navigate (or used to be). Both have the same data, but are better for the extras than NCBI.


Have DNA motifs 6-12bp long, trying to get conservation scores

In case you have only ACGT in your motifs

The short motifs make it sound as if you are in the business of looking for a kmer counter. You can choose to use existing software or build your own.

Using existing software might be your easiest path. An older post from 2014 will probably give you a first idea what's out there: http://homolog.us/blogs/blog/2014/04/07/kmer-counting-a-2014-recap/ . Note that a couple of algorithms mentioned there got successors, so it is worthwhile digging a bit around. The small kmer size will make most of them usable for your needs.
As the maximum size of your kmers is comparatively small (12 nt need 24 bits, i.e., max 16.7 million entries in your kmer table), you should be able to easily roll your own kmer counting in about any language you like and on about any of nowadays computer. The pseudocode section on the Wikipedia entry for kmers will give you first pointers for that. Might be a bit more work, but maybe more flexible depending on your needs.

In case you have IUPAC bases (N, W, etc.) in your motifs

I don't know any any pre-existing software doing what you need. I could imagine that the short motifs make using regular expressions doable for this kind of search, but I may be wrong. Testing this should be easy in a simple script as all major programming languages have modules or libraries for REs. Even if it should take a couple of hours to run on your data set, that would be good enough for a one-off calculation.


What is the detection limit of protocol in scRNA-seq

The detection limit refers to the ability to detect a given gene/transcript. In this case, they're saying that a given transcript requires 5-10 copies in order to be likely detected (in other words, the loss rate is 80-90%, which is directly related to a previous question of yours).


how to have a blast xml file in a readable and understandable way

The most human-readable blast output formats are 0-4, e.g. -outfmt 0 - which if I am not mistaken is the format used by web blast. However, these formats are a pain to automatically parse.
As Cowboy_Patrick pointed out, XML is more of a machine-readable format. The tabular output format ( -outfmt 6) is very commonly used, because it is easily parsed by unix command-line tools, and its structure also allows a somewhat easy human reading and understanding.


You seem to refer to the GSEA provided by the Broad institute, (there are other GSEA algorithms).

1) You can provide whatever you wish, but if you want to know if those gene sets in which side of the ordered list are they, then provide all the list (of genes) you have.

2) GSEA analyize if the order of a given list distributes in certain way the elements of another (unordered) list. In your case you would need to compute the GSEA for each treatment vs control comparison.


COG Annotation - Dealing with genes assigned to two or more COG categories

You add 1 to the count of each COG category. You are looking for over represented COG categories so you must count them all. Many genes will be assigned to multiple categories. In fact, a majority of human genes are given multiple GO cellular function annotations. I have worked considerably in this field and one of the surprises I found is how common protein multifunctionality is. But I digress.

The point here is that you want to see whether any COG categories are over represented in your test set. Whether a gene is annotated to one or to a hundred categories is irrelevant. You are counting genes per COG, so you need to count each of the COGs assigned to a gene separately.


determine if ChIP-seq peaks are broad or narrow

I'd say mostly this is a question of understanding the underlying biology and the relevant literature.

If it is not known in the literature whether a mark is peaky or broad, evidence might come from FISH studies or Low throughput qPCR. Another way to look would be to examine the signal expressed as fold enrichment over input on a genome browser and look if peaks tended to be small (a few hundred bp, e.g. a transcription factor), or large (several 1kb, e.g. H3K27me3) or contain both (e.g. RNA polymerase II).


How can I use annotations to remove variants not relevant to cancer risk?

While your question is specific to cancerous germline mutations, I'd suggest you look at the COSMIC database of somatic mutations to include in your analysis.

There are other factors to include in this kind of analysis you're suggesting, such as predictive deleterious effects (PolyPhen for example can perform such predictions).

If you have 10M variants/indels, for the same form of cancer, then look for common variants, or maybe the frequency of variants identified across exons of a gene.

What you've asked is highly broad and we'd need more specific details about your cohort to suggest more specific approaches.


What is the best way to drop alignment artifacts in sequence objects in biopython?

I think the .ungap() is what I was looking for.
def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.ungap().translate())} for rec in fasta_sequences])
    return sequences


You need to include bwa.o in your call to clang, since it holds the definitions of those functions. This will cause a good number of object files, so you might as well just clang -o foo *.o foo.c -lz -lpthread -lm -lrt and get it over with. Since this was quite obviously never meant to be done, you'll have to remove at least one main() function from bwa.c.

Note that further questions regarding this will be considered off-topic and closed, since it's purely about compiling something that you want to write yourself.


validating identified sub-populations of cells in scRNA-seq

A SC3, single-cell consensus clustering, approach could be helpful here. It aims at achieving "high accuracy and robustness by combining multiple clustering solutions through a consensus approach" https://www.nature.com/nmeth/journal/v14/n5/full/nmeth.4236.html


Answer from @throckmorton, converted from comment:
Here is a recent approach to this problem, using multiple linear references:
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02229-3


How are positions visually indicated in the UCSC Genome Browser?

It's an A. The number refers to the position just before the |. This is clearer on different zoom levels, but I double checked by going to the very beginning of a chromosome:

As you can see in the image above, the | is always after the base it refers to. I also confirmed by checking a few random positions using another genome browser which places the positions unambiguously above each base.

For more details, see this UCSC blog post.


How should the popular press compare similarity of genomes?

1–4% is from an evolution point of view. 99.7% is from a sequence point of view. Because they are measuring different things, they can be compatible with each other. The correct interpretation is: 1–4% of a non-African genome is inherited from Neanderthals and the sequence of this 1–4% differs from modern human sequence by 0.3%. PS: so, for each non-African, Neanderthals contributed 1–4% * 0.3% of base pairs not present in modern human. This sounds tiny, but even a single SNP may lead to severe disease. This tiny fraction has known functional impact to diabetes and high-altitude adaption.

Others pointed out two related questions:

How many genes do we share with our mother?
Do apes and humans share 99% of DNA or 99% of genes? What is the difference?

These questions are really about the ambiguity in the definition of "gene". Importantly, traditionally in popular science and population genetics, a gene could mean an allele. This makes differences.

For example, by common uses of "gene" in biology, the answer to "how many genes do we share with our mother" would be nearly 100% – every human has nearly the same set of genes. However, clearly this is not what the asker wants. By gene, he/she really means allele. The expected answer is 50% – we inherit 50% of alleles from our mother.

In the end, what is a "gene" exactly? I don't know. Biologists can often figure out the meaning of word "gene" in the right context, but non-biologists will sometimes have troubles with the ambiguity.


Get protein names corresponding to PDB ID

You can use one of the UniProt Protein APIs.
As you said you have your pdb entries in a text file line by line you can, like this example.txt containing:
1brr
4lzm
2dyi

Using the commandline, you can use a little script like this to download the name, if it is available for the given pdb entry.
while read line;
do
 curl -X GET --header 'Accept:application/json' "https://www.ebi.ac.uk/proteins/api/proteins/pdb:$line" |
 jq -r '.[].protein.recommendedName.fullName.value'  |
 sed "s/^/$line\t/" >> pdb_names.txt;
done < example.txt;

You need to have curl, sed and jq installed on your system.
This gives you following output in pdb_names.txt
1brr    Bacteriorhodopsin
4lzm    Endolysin
2dyi    Ribosome maturation factor RimM

Update
if you want to speed it up, you can run it with parallel.
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" |   jq -r ".[]. .protein.recommendedName.fullName.value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

With the -j option you call how many jobs should run in parallel. The limit of the UniProt API is 200 request per second per user.

Update 7. Nov 2020
To get another info beside the protein name, you need to know how the JSON response from UniProt looks like.
To get also the scientific name, you can run following command:
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" | jq -r ".[] | .protein.recommendedName.fullName.value + \" - \" + .organism.names[0].value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

As result you get this:
1brr    Bacteriorhodopsin - Halobacterium salinarum (strain ATCC 700922 / JCM 11081 / NRC-1)
4lzm    Endolysin - Enterobacteria phage T4
2dyi    Ribosome maturation factor RimM - Thermus thermophilus (strain HB8 / ATCC 27634 / DSM 579)


How to fit all genes (labels) in chromosome ideogram plot made by RCircos package?

To change the maximum number of genes in RCircos, we could modified the char.width according to this link.
The actual value in a RCircos session could be modified with get and reset methods for plot parameters.
params <- RCircos.Get.Plot.Parameters() #$char.width
params$char.width <- 100 #default 500
RCircos.Reset.Plot.Parameters(params)

#the maxLabels are updated accordingly
RCircos.Get.Gene.Name.Plot.Parameters()


Why can using more reads lead to a lower quality assembly?

The more reads you add the more errors you add into the assembly. This is because additional duplicate reads don't add nodes/edges to the de Bruijn graph, but those with errors do. By preselecting those reads aligning to your desired contigs, you may be further exacerbating this effect and further hindering SOAPdenovo2's ability to do k-mer correction (after all, the observed k-mer frequency of errors would be higher). This then results in additional incorrect paths through the graph that the assembler has to traverse.


How to extract all peptide sequences for a UniProt protein entry with an API?

If downloading the complete database is not an option, you may able to use the instructions from this help page (based on the gff format for any given query):
https://www.uniprot.org/help/download_subsequences


How to release an R package with several bed files?

You can store it in data/ subfolder. If one of your functions need this bed file, you can include importing it in your function.

Here are some examples how to get the data from data/ subfolder.

There are probably more ways to do this, this is just one way.


What kind of BLAST do I need to do to accomplish this task?

Clarification
Just making sure - you have ~6000 datasets - one for each S.cerevisiae gene of interest - which are made up of homologues of each respective gene? And you want to filter out genes that are too similar?
In that case, it would be helpful to know how you're defining a sequence as 'similar enough' for it to be filtered. Given they're all homologues, it's likely there's going to be a high baseline of sequence similarity. In your example, you consider 97.40% identity across the whole query as different enough to be filtered out, so I'm getting the feeling you specifically want to filter out identical sequences?
Also, when you talk about genes here, do you want to filter out gene products (i.e. proteins) that are too similar to the actual coding sequences (CDS)? The CDS  of two homologues might be X% different, but that difference might not be as high once translated into proteins due to how codons work. Proteins can go down to levels of similarity as low as ~30% identity and still have the same function, so depending on the context of your work (taxonomic, structural, functional) you might want to translate these CDS into protein, so further clarification would be helpful!
The actual bioinformatics question
To answer your question of programming, if the context of this comparison is always X S.cerevisiae gene vs. Y dataset of homologues, a MSA is definitely not the answer (you would do that if you wanted to see how all of the homologues compare to each other, rather than to specifically X gene). Also, using the web-client for NCBI BLAST will take a lot of time, so that's also not recommended for the kind of throughput you need.
Pairwise sequence alignments that you do on your local machine/server between X genes and all of its respective Y dataset sequences is probably the best way to go. I would suggest looking into needleall - it allows you to compare X sequences to Y sequences in a pairwise way, though it performs a global alignment rather than local like BLAST, so it might be slower.
I don't know what your computational resources are, so in case that might be too computationally heavy, the faster alternative would be to install the BLAST+ suite, which allows you to run BLAST queries like you did on the web browser, except on your own machine. Specifically, you could make your own custom individual BLAST databases out of each Y dataset using makeblastdb from the suite. Then you use blastn (or blastp if you decide to translate your sequences) on each dataset with its respective S.cerevisae sequence.
BLAST+ also allows you to customise the actual output as well, making it easier to parse your results programatically for when you want to remove sequences that are too similar. There are many examples online of how to do this, but have a go at it on your own  first with the BLAST+ documentation I linked :)


R package development: How does one automatically install Bioconductor packages upon package installation?

As suggested, here’s an example showing the relevant lines from a DESCRIPTION file from a CRAN/GitHub hosted project that has Bioconductor dependencies (truncated):

Depends:
    R (>= 3.3.0)
biocViews:
Imports:
    methods,
    snpStats,
    dplyr

The relevant bit is the empty biocViews: declaration, which allows the Bioconductor dependency {snpStats} to be automatically installed.


Looking for genes that are expressed the same in all the samples

I cannot think of a way without performing differential expression analysis. Since you don't have controls, what you could do is assigning some of the samples as controls.
Let's say you have samples A, B, C, D, E and F. You can perform differential expression analysis A, B, C (i.e. treatment) versus D, E (control). If you perform this for all possible combinations (that would be 5 choose 2 which is equal to 5 choose 3 if you decide to go for 3 vs 2 splits) and filter for genes which do not show up in any DE list you end up with the genes that are expressed similarly in al of these samples.


Yes, if you load the list you are interested in as a non-indexed "BED" file then you can jump to an rsID or any other name in the BED file. You'll need a targeted list, dbSNP for example is much too large to load in this way.


Ensembl protein identifiers from different assemblies

Those IDs are elderly! Ensembl 54 was 2009!

I would recommend using BioMart combined with the ID history converter. The ID history converter will convert old IDs to new, and BioMart will convert ENSPs to ENSGs. You can either use the ID history converter with the ENSPs first, then the current BioMart to get the ENSGs. Or you can use Ensembl 54 BioMart to convert the old ENSPs to old ENSGs, then update the ENSGs with the ID history converter. I would probably do both and compare results.


Thanks for pointing this out. Looks like we forgot to include two important lines. Something along the lines of

ConnectivityHuman = colSums(adjHuman)-1
ConnectivityChimp = colSums(adjChimp)-1

This should do it; I'll update the tutorial.


FPKM are inherently experiment specific and can not be used to compare across samples. Let's consider the following two sequencing runs. Let $E1$ and $E2$ be the true, underlying expression in two samples of genes 1-6. Let $S1$ and $S2$ be the observed expression in our sequencing.
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 10 & 100 & 20 \\
  G2 & 100 & 10 & 100 & 20 \\
  G3 & 100 & 10 & 100 & 20 \\
  G4 & 100 & 10 & 0 & 0 \\
  G5 & 100 & 10 & 0 & 0 \\
  G6 & 100 & 10 & 0 & 0 \\
\end{matrix}
$$
Our totals are 60 counts for $S1$ and 60 for $S2$. We sequenced both libraries to the same depth, but as there are fewer genes expressed in $E2$, we just naturally capture the genes that are expressed more often. For simplicity's sake, let's assume the genes all have the same length of 1, so we ignore the "K" in FPKM. Let's also forget the scaling by a million, since that's just to make nice numbers:
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 0.167 & 100 & 0.333 \\
  G2 & 100 & 0.167 & 100 & 0.333 \\
  G3 & 100 & 0.167 & 100 & 0.333 \\
  G4 & 100 & 0.167 & 0 & 0 \\
  G5 & 100 & 0.167 & 0 & 0 \\
  G6 & 100 & 0.167 & 0 & 0 \\
\end{matrix}
$$
Now let's calculate size factors. In DESeq2, we create a pseudo-sample which is the geometric mean of all counts of each gene, so $\sqrt{10 \cdot 20} = 14.14$ for genes 1-3 and genes 3-6 are ignored as they have a member with 0-valued counts. The size factor for each library is the median ratio for each library compared to that pseudo-sample:
$$
med([\frac{10}{14.14}, \frac{10}{14.14}, \frac{10}{14.14}) = 0.71
$$
$$
med([\frac{20}{14.14}, \frac{20}{14.14}, \frac{20}{14.14}]) = 1.41
$$
So where do we end up if we divide by the size factors?
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 141 & 100 & 142 \\
  G2 & 100 & 141 & 100 & 142 \\
  G3 & 100 & 141 & 100 & 142 \\
  G4 & 100 & 141 & 0 & 0 \\
  G5 & 100 & 141 & 0 & 0 \\
  G6 & 100 & 141 & 0 & 0 \\
\end{matrix}
$$
So in this simple example, size factors produce a better estimate than FPKM or TPM to compare values across experiments, hence the "robustness". Please note that this is still not great for any statistical testing, and in general should be used with caution anywhere.
Finally, I'll leave some more reading material concerning FPKMs and cross-sample comparisons here:

[1] Dillies M.A., Briefings in Bioinf. 2012
[2] Soneson & Delorenzi, BMC Bioinf. 2013
[3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: https://www.youtube.com/watch?v=5NiFibnbE8o&t=32m
[4] From Answer: "Over-correction" in the size-factors of the DESeq2 package: https://support.bioconductor.org/p/80733/#80736
[5] http://www.biomedcentral.com/1471-2105/14/370
[6] https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/
[7] http://rnajournal.cshlp.org/content/22/6/839


Merge multiple PDB files into a single file

PyMOL is a rather flexible tool for doing this.

The command create combo, domA or domB or domC in the PyMOL app will create a single object. However, they will need to be the same chain and segment (segi). Also you will need to fix the numbering to match what you expect them to be.
To fix this you need to do:

 alter domA, chain='A'
 alter domA, resv+=0
 alter domB, chain='A'
 alter domA, resv+=100
 alter domC, chain='A'
 alter domA, resv+=200
 alter all, segi=''
 sort
 create combo, domA or domB or domC
 save combo.pdb, combo

PyMOL can be installed as a Python 3 module, which works in the same way but with pymol.cmd and is a lot more replicable. But for starting out the App is fine.

import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.alter('domA','chain="A"')
     etc.

Now. This assumed that the coordinates are where we want them to be. Which is highly unlikely and that there are zero overlaps.

If there are overlaps, aligning the overlapping residues (align domB and resi 80-100, domA and resi 80-100) can be done. Just remember that the first element is the mobile object, the second the target, hence why I wrote it backwards.
If do not have overlaps and want the domains docked together you can use the DEMO tool by the I-Tasser people.
If you want your domains as "balls on a string" aligned on the x-axis for a nice picture and you do not care if it may not be biological, you may want to check out a Github repo of mine for doing just so (it offsets, aligns if possible and spaces them out, but does not add linkers).

If you want to add linkers, you can use submit your sequence to Phyre or I-Tasser with your missing density model as a reference, or you can add them with rosetta remodel (tricky) or make them up in Coot or PyMOL (bad idea, but if you are in a hurry... here is a tutorial of mine).


how can i search in p-blast against UNIPROT(not only swissprot)

I don't know of any online tool enabling a search against TrEMBL. You probably will have to set up your own local BLAST environment, roughly following these steps:

download blast+ - but there might be an easier way to install the software, depending on your operating system (e.g. on Ubuntu you can just sudo apt install ncbi-blast+)
download target databases
build local database

I hope this helps.


Detecting differentially expressed genes with foldchange >= 2 and FDR < 0.05

The edgeR authors recommend that you use a relatively low logFC threshold for glmTreat such as lfc=log2(1.2). A lfc value as high as lfc=log2(2) is seldom required and only would only be appropriate for datasets with really large amounts of DE.
The first threshold you used of lfc=2 is equivalent to lfc=log2(4) and is too high for any dataset.
There's a few subtleties here:

First, there's the difference between logFC and FC as explained by Devon.

The lfc in glmTreat is not a fold-change cutoff! It is instead a testing threshold. All the significant genes will have actual logFCs well above the threshold.

The lfc threshold applies to the raw fold-changes (unshrunk.logFC in the output), not to the shrunk versions (logFC in the table).

The intent of glmTreat is to give a better gene ranking than would be provided by FDR or fold-change, so please don't re-sort the results by fold-change -- that would tend to defeat the purpose of the method.

I understand the documentation for glmTreat is a bit light, but you can get an idea of how to use it from our workflows.

https://bioconductor.org/packages/release/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html

https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html

I understand the simple fold-change cutoff of FC > 2 with FDR < 0.05 is often used in the literature, but the edgeR authors do not recommend it and hence we don't provide an automatic way to get it from the edgeR pipeline.
For RNA-seq, fold-change cutoffs tend to prioritize genes with low expression and low counts.
The glmTreat approach is more subtle but we think gives better results, both in terms of biological relevance and statistical correctness.


What is the difference between a transcriptome and a genome?

In brief, the  “genome”  is the collection of all  DNA  present  in  the  nucleus  and  the  mitochondria of a  somatic  cell. The initial product of genome expression is the “transcriptome”, a collection of RNA molecules derived from those genes.


How to trigger next process on completion of previous process?

So from this page...
https://groups.google.com/g/nextflow/c/jt77_-uApMs
It is the "scatter-gather" pattern. I think the answer is to collect() the channel as follows - interested to hear if there is a better way
files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    val go from result.collect()

    script:
    """
    echo "Triggered once after all files complete!"
    """

}


Error is the estimation that can't be fitted by your parameters (in this case g).
The amova function in pegas package does the following:

  This function performs a hierarchical analysis of molecular variance as described in Excoffier et al.
  (1992). This implementation accepts any number of hierarchical levels

And as per the definition of Fixation index, you need the variance in the frequency of the allele between different subpopulations, weighted by the sizes of the subpopulations, the variance of the allelic state in the total population. Which is not provided by amova. Thus you cannot compute the fixation index from amova results.


40% missing data is huge. Missing data analysis is complicated on the underlying distribution. If the data set is periodic then missing data periodicity is needed. Non-periodic data can be solved using regression, I think SciPy has an automated method for this, but the volume of missing data makes this approach complex. If you can bring the missing data down to 20% by compartmenatlising the data set that would be the starting point.


Is it possible to obtain RefSeq transcripts from Ensembl database for GRCh37 to match Ensembl Genes?

The RefSeq match option in BioMart is from the Matched Annotation from NCBI and EBI (MANE) collaboration between RefSeq and Ensembl. It has only been calculated for the up-to-date gene annotation on GRCh38 so cannot be obtained on GRCh37. You can get mapping from Ensembl to RefSeq transcripts through BioMart as RefSeq mRNA ID (refseq_mrna in R) but this is not a perfect match like the MANE, it is a mapping based on sequence similarity and similar genomic location, and there can be mismatches between them.


UCSC liftover cannot produce bed format output

I figured out the issue by myself. It turned out that UCSC liftover doesn't like the character ":"

I removed ":" and it worked


snakemake with --conda-prefix, can I use a pre-built conda environment?

If the environment already exists it will simply be used. You cannot, however, tell it to use an environment with a normal name (e.g., "my_env").
I assume you mean, "where the heck does the random looking conda environment name come from?!?", which is a very good question. The answer to that is that it's a hash of the --conda-prefix setting and the contents of the environment yaml file you specify. The actual code for doing this (I'm "stealing" it from snakePipes) is:

import hashlib
md5hash = hashlib.md5()
md5hash.update("what you gave to --conda-prefix".encode())
f = open("your environment yaml file", 'rb')
md5hash.update(f.read())
f.close()
h = md5hash.hexdigest()

h is then the directory name. There are actually two lengths of those, the shorter 8 character version that snakeMake uses by default or the full length version that the code above will produce. snakeMake will actually use either of these. If you're wondering, no I don't think any of this is documented anywhere, you have to dig through the snakeMake code to find it.


In the dbSNP FTP site you linked, you need to go into organisms and select your organism of interest (human obviously). You can then select which release of the SNP database you want (eg. 150 vs 151) [the modified dates hint at when they were released] and which genome build you want (eg. hg37 vs hg38).

Below I am looking at the .vcf files. I assume this what you want because they can be converted to ADAM format via the adam vcf2adam command.

For example human (newest genome and newest dbSNP release):
ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/

Your next decision is whether or not you want to use ALL the annotated SNPs in the database or just the common ones. For the former you choose the All_20180418 file and for the later the common_all_20180418 file. I'm pretty sure the convention here is that the date in the filename is when it was uploaded/updated.

In your case I'm not sure what makes the most sense. Firstly, are you trying to replicate the initial conditions of when that software was published or do you want simply to benchmark the method with the best available data. Besides that, I imagine using all them is computationally much more taxing and if you are just coming up with general trend then using the common ones seems viable.


RNA velocity: competing explanations for variable ratio of spliced to unspliced transcript

If cells aren't in a steady state then one wouldn't expect to see the results in (e), where there are two slopes. One would instead expect to see results more like Plekha3 in panel (c), where there's more variability around the the steady state line (or perhaps no steady state line at all).


WGCNA correcting p value for varying sample numbers in each of the trait columns

You can use the function corAndPvalue. Simply replace the lines

moduleTraitCor= cor(MEs, status, use="p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

by something like

cp = corAndPvalue(MEs, status)
moduleTraitCor = cp$cor
moduleTraitPvalue = cp$p

corAndPvalue takes automatically into account the number of valid samples for each element of the resulting correlation matrix.


Coverage calculation: long reads (RNA-seq)

The read length is irrelevant when calculating the mean coverage statistic. It's simply the total number of bases sequenced divided by the target Xome length.

In the example provided in the question, $L*N$ is equivalently expressed as $\sum^{N} l_x$, or the total length of sequenced reads. For equal-length reads, this is easier to calculate as number of reads $*$ length of reads, but the full sum is necessary if read lengths are not equal.

For what it's worth, this value, $L*N$, is also equivalent to the average read length multiplied by the number of reads, i.e.

$$\frac{\sum^{N} l_x}{N} * N$$

... which makes it a bit more obvious why the total length of sequenced bases is very slightly easier to calculate.

In answer to the second part of the question, $N$ should ideally correspond to only the sequences that are from the reference. For an assembly, only the reads that were used to create the assembly should count when considering coverage. When there's no likelihood of contamination (and especially for RNASeq with variable coverage per transcript), a statistic using the total number of reads should be good enough for most purposes that I can think of. High-throughput sequencing is such a fuzzy process that a few tens of millions of bases here or there shouldn't make a difference.


Bowtie2 is no longer the fastest aligner. Salmon and Kallisto are much faster, but have been designed to optimise RNASeq mapping. Their speed is gained from avoiding a strict base-to-base alignment, but they can output mostly-aligned reads (i.e. position-only, without local alignment) as pseudo-alignments. See here for more details.
Both Kallisto and Salmon can do additional bootstrapping that is interpreted by sleuth for improved performance in isoform detection. They can also output counts that are equivalent to read-level counts from other programs, which can then be used by other downstream gene-based differential expression analysis software (e.g. DESeq2). Salmon has additional options that can  correct mappings for sequence-level and GC bias.
HISAT2 is from the same group as Bowtie2, and does the same sort of stuff, but with a few optimisations added on top. In particular, it's much better at working out split reads from RNASeq runs. Like Bowtie2, it will do local alignment of reads.
For quick alignment of long reads, minimap2 works well. For high-accuracy alignment (but comparatively slower), LAST works well.
Most bioinformaticians seem to prefer STAR for things that Bowtie2 was previously used for. I'm not yet convinced it's a better alternative, and currently prefer HISAT2 for high accuracy short-read alignment.


Instead of outputting all files and filtering for only bim, output only bims:
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)

Do
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into step1a_out
tuple val("${prefix}"), path("${prefix}.bim") into step1a_bim

So you can skip your complicated filter:
process step1b {
    input:
    tuple val(bim), path(bimfile) from step1a_bim
    output:
    path("${prefix}.step1b.snplist.txt")
tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into
step1b_results
    script:
    """
    awk '{ if ((\\\\\$5=="T" && \\\\\$6=="A")||(\\\\\$5=="A" && \\\\\$6=="T")||(\\\\\$5=="C" && \\\\\$6=="G")||(\\\\\$5=="G" && \\\\\$6=="C")) print \\\\\$2, "ambig" ; else print \\\\\$2 ;}' "$bimfile" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"

"""

```


The counts files for GSE89225 is the output of HTSeq-count as a large matrix. Unless you are developing a differential expression package yourself you should not attempt to directly use this. Rather, you should load it into R and use packages such as DESeq2, edgeR, or limma (those are the most popular ones).

For convenience, in DESeq2 you would want the DESeqDatasetFromMatrix() function after loading this.


I had a similar case just recently.
The tssRegion argument ChIPseeker is used to define what is considered a promoter, which is frequently the question at hand. However, ChIPseeker will not perform any kind of filtering on the annotated peak file.
In the results there is a column annotation that contains information what kind of genomic region the peak is in. If it is within your tssRegion window it will be labeled as a promoter.
In your case the annotation is 1st Exon = Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9. This tells you that the peak is not part of the promoter region you defined and located within the first exon of the indicated transcript. The peak is really far away from TLR4 (374699 bp) but I think all other sequences nearby are pseudogenes and this is why TLR4 is selected as the closest gene.
You can get an overview over the annotations using plotAnnoPie(peakAnno)
tl;dr: you need to filter the results of the ChIPseeker annotation to get the genomic regions you want.


Output the linear predictor from a stratified cox model?

Unlike wikipedia, the linear predictor (LP) in coxph is made dependent from the baseline hazard function $\lambda$. For example, if a model has two predictors x1, x2 stratified by sex, the linear predictor would be like:
$$LP = log (  \lambda _{sex}) + \beta _{1}x _{1}+\beta _{2}x _{2}$$
Where the $log(\lambda _{sex})$ can be obtained by predict() with x1=x2=0.
We could calculate the LP without strata information (sex) by hand using the $\beta _{1} \, \beta _{2}$ provided by coef(fit) and a customised $\lambda _{unkown \, sex}$.
Below is a simulation comparing the results from predict() and by hand coef(fit).
## Toy data
df <- list(time=c(4,3,1,1,2,2,3),
          status=c(1,1,1,0,1,1,0),
          x1=c(0,2,1,1,1,0,0),
          x2=1:7,
          sex=c(0,0,0,0,1,1,1))

## Fit a stratified model
fit <- coxph(Surv(time, status) ~ x1 + x2 + strata(sex), df)
coef(fit)
#         x1          x2
# 0.79451881 -0.01917633

## Baseline linear predictor
lambda_sex1 <- predict(fit, newdata=list(x1=0, x2=0, sex=0))
# -0.746578
lambda_sex2 <- predict(fit, newdata=list(x1=0, x2=0, sex=1))
# -0.1497816

## by predict function
predict(fit, newdata=list(x1=1, x2=2, sex=0))
# 0.009588167
predict(fit, newdata=list(x1=1, x2=2, sex=1))
# 0.6063845

## by hand
lambda_sex1 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.009588167
lambda_sex2 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.6063845

Finally, to calculate "a" linear predictor for an unknown sex:
# prevented by predict function
predict(fit, newdata=list(x1=1, x2=2, sex=2))
#Error in model.frame.default(data = list(x1 = 1, x2 = 2, sex = 2), formula = ~x1 +  :
# factor strata(sex) has new level sex=2

# by hand, assuming lambda for sex 2 is 0
0 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.7561661


How to download SNP data from specific regions, population and positions?

The HapMap website was suspended due to a security issue. Archived HapMap data is available via FTP, but it's a better idea to download up-to-date International Genome data from the IGSR website (which includes the HapMap individuals as a subset).

The format you use will be dependent on the data analysis program you're using. In most cases involving SNP analysis, the VCF format is probably what would be used.


The question of folding of RNA sequences looks slightly similar to protein folding - perhaps searching in this domain might bring more suggestions.

An example of (current) state of the art of deep learning approach to protein folding is AlphaFold developed by DeepMind for CASP competition . They basically decompose the angles to an image representation and look for structures in the proximity of each atom. While it requires a significant processing power, RNA structures might be easier due to the size of the sequence.


Should the cell sorting marker genes be excluded during clustering?

I do not think there is a simple "yes" or "no" answer here.

A good starting point would be, as you suggest, use all the genes and assess the results in the light of the marker genes and expected results. This could both serve as as good quality control as well as give you overview of all the processes happening in the cells.

Depending on the marker genes effect and biological question you may then want to remove the marker genes with potential ordering effect, or even other genes, e.g. by GO terms or pathways.

You most likely want to account for the cell cycle phases as well. And check for other technical factors. I recommend Bioconductor pipeline to get inspired https://www.bioconductor.org/help/workflows/simpleSingleCell/ when it comes to scRNA-seq analyses.


You could try aligning your reads to the draft reference genome with for example minimap2 and calling variants with freebayes. It appears there is a protocol for long reads: https://github.com/ekg/freebayes#re-genotyping-known-variants-and-calling-long-haplotypes

Disclaimer: I have not run freebayes myself, and it may be that this is a terrible idea. I would expect that calling SNVs (Single-nucleotide variants) from PacBio reads is very challenging.  However, with enough coverage, I think it could be possible.


>My counter is counting genotypic combination occurences more than once, how do I ensure it counts one combination and doesnt count it again?

Alright, so there are a number of problematic patterns in your code - as far as I understand what you are trying to do. Next time, try to post a reproducible example that people can use and more people will be willing to help.
L18-21:
combination_labels = []
combination_counts = []

for lineage in lineages:

Declaring these two lists before the loop, then appending all your combination labels inside the loop means that these two lists will contain all labels and counts for all runs (all your lineages). If I am understanding your intentions correctly this is probably not something you want. The simple fix here is to just move them into the for-loop.
L29-30:
for mutation in combination:
    combination_df = combination_df[combination_df[mutation] == 1]

Here you are iteratively checking if the value for each mutation is set to '1', but you are also overwriting your combination_df variable at each iteration of the loop. After this loop, combination_df will be whatever mutation last was set to '1' in your current combination. I will get back to solving this in my solution below.
L37-44:
combination_labels.append('_'.join(combination))
combination_counts.append(combination_df.shape[0])

out_df = pd.DataFrame({'combination':combination_labels,
                   'count':combination_counts})
out_df['percentage'] = (out_df['count'] / df.shape[0]) * 100
out_df = out_df.sort_values('percentage', ascending = False)

out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv',
         header = True,
         index = False)

Since this part is again scoped inside the for-loop, each iteration (each combination) will overwrite the same 'Results_2/' + lineage.replace(".csv", "") + '_3.csv' file. This you also want to move out of the for combination in ... loop.

A big concern is also performance. You have 2000-something combinations and 100,000 lines in your file. As you loop over combinations, you also process the DataFrame each time, when this is a proble you can solve in one iteration. It's not much data so it's still fine, but these habits are still good to develop early.
Now here's an idea for a rewrite that should work better. I don't have a csv from you to test with, so there might (proabably are) still things to fix, but hopefully it'll give you a starting point:
mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

combinations = []
for M in range(1, len(mutations)+1):
    for subset in itertools.combinations(mutations, M):
        combinations.append(subset)

combinations = ['_'.join(sorted(x)) for x in combinations]

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

for lineage in lineages:

    out_df = pd.DataFrame({"Count": 0}, index=combinations)

    df = pd.read_csv('Results/' + lineage, index_col="lineage").dropna()

    for index, row in df.iterrows():
        combination = '_'.join(sorted(df.columns[row==1]))
        out_df.loc[combination, "Count"] += 1

    out_df['Percentage'] = (out_df['Count'] / df.shape[0]) * 100
    out_df = out_df.sort_values('Percentage', ascending = False)

    out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv', header = True,
        index = True)


What are the differences between BNGL (BioNetGen Language) and SBML (System Biology Markup Langague) formats?

Both formats have formal grammars for describing similar concepts, however they are different and to my knowledge one is not a proper subset of the other. While one could do a comparison of a syntactic nature, a simple Google search of just both acronyms reveals BioNetGen 2.2: Advances in Rule-Based Modeling

In here can be found:

  2.4 SBML-to-BNGL translation
  SBML is a widely-used model exchange format in systems biology (Hucka et al., 2003). Models encoded in SBML
  are ﬂat, i.e., their species do not have internal structure, which
  limits their utility for rule-based modeling. BioNetGen 2.2 includes
  an SBML-to-BNGL translator, called Atomizer (also available as a web
  tool at ratomizer.appspot.com), that can extract implicit molecular
  structure from ﬂat species (Tapia and Faeder, 2013). A full report on
  Atomizer and its application to the BioModels database (Li et al.,
  2010) is currently in preparation. However, Tapia and Faeder (2013)
  reported that an early version of the tool could recover implicit
  structure for about 60% of species in models within the database that
  contain ≥20 species. Thus, Atomizer makes available a large set of
  pre-existing models in a rule-based format, facilitating their
  visualization (Wenskovitch et al., 2014) and extension (Chylek et al.,
  2015).

While not directly related to the question this also helped in understanding more about BioNetGen

Wikipedia: Multi-state modeling of biomolecules

  Biological signaling systems often rely on complexes of biological
  macromolecules that can undergo several functionally significant
  modifications that are mutually compatible. Thus, they can exist in a
  very large number of functionally different states. Modeling such
  multi-state systems poses two problems: The problem of how to describe
  and specify a multi-state system (the "specification problem") and the
  problem of how to use a computer to simulate the progress of the
  system over time (the "computation problem"). To address the
  specification problem, modelers have in recent years moved away from
  explicit specification of all possible states, and towards rule-based
  formalisms that allow for implicit model specification, including the
  κ-calculus, BioNetGen, the Allosteric Network Compiler and others. To
  tackle the computation problem, they have turned to particle-based
  methods that have in many cases proved more computationally efficient
  than population-based methods based on ordinary differential
  equations, partial differential equations, or the Gillespie stochastic
  simulation algorithm. Given current computing technology,
  particle-based methods are sometimes the only possible option.
  Particle-based simulators further fall into two categories:
  Non-spatial simulators such as StochSim, DYNSTOC, RuleMonkey, and
  NFSim and spatial simulators, including Meredys, SRSim and MCell.
  Modelers can thus choose from a variety of tools; the best choice
  depending on the particular problem. Development of faster and more
  powerful methods is ongoing, promising the ability to simulate ever
  more complex signaling processes in the future.


I had a think about this one and am stumped. I thought there could be a connection with acute oak disease (AOD) or else chronic oak disease (COD), which are emergent infections threatening the population. However, they appear to be fungal and bacterial infections rather than viral.
Single stranded RNA plant viruses are very common, more so than in mammals, but to achieve the results you are obtaining the trees would need to be really sick if the predominent RNA was viral. Short viral RNAs are produced regularly by viruses as defective RNA and its a major area of research, considering a mechanism of defective interference,
In context it would be weird if this RNA swamped the tree/plant, you need an actively replicating virus for this stuff to appear. So I don't think its either viral or defective viral RNA.
The antiviral response in mammals can involve RNase and this could explain your degradation. Essentially, mammals (probably vertebrates) kick-back hard when virus appears. However, I couldn't really find evidence of this in plants.
I don't think the results you obtained have a clear biological explanation in pathogenesis. The best explanation is it was wet lab-technique, e.g. library construction resulted in overdigestion/fragmentation. The only thing I can suggest is template assembly using a normally distributed plant/oak RNA library.
I was probably overthinking in forging a link with an emergent (tree) pathogen.

!


How do I re-name the headers of my Fasta file?

If awk is good as well:

$ awk -v FS="|" -v OFS="|" '/>/{
    gsub(";",",", $5);
    gsub("__",":",$5);
    print ">"$2,$3";tax="$5"_"$3";";
    next;
}1' file.fa
>JF514531|SH174817.07FU;tax=k:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

This relies on the order of the parts separated by |.

I don't understand why k__Fungi becomes d:Fungi, so I ignored it for now.

Some explanation for what awkis doing.

| is defined as a field (column) separator for the input and output
If the line contains a > the code between { and } is executed.
the first gsub replaces ; with , in the 5th field
the second gsub replaces __ with : in the 5th field
then we print the new header line. The $2, $3 and $5 get replaced by the values of the second, third and 5th field of the input. , get replaced by the defined output field separator
with next we tell awk to go to the next input line and repeat the code
The 1 is in equivalent to print $0 for printing the whole input line. In our code this is only executed if the line is no header line.


Challenging benchmarks for supervised learning on sparse scRNA-seq data

You could try the PBMC 3k data from the Satija lab:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
For well-annotated data, there's the Single-cell proteo-genomic reference map:
https://www.biorxiv.org/content/10.1101/2021.03.18.435922v1
Data for that are available here.
If these datasets are too large (i.e. too good), you can use the existing counts as weighting for a random removal of reads. Genes with high expression are likely to be present regardless of sequencing depth, hence the count-adjusted removal.
FWIW, I don't think that supervised learning is a good fit for dealing with data sparsity in cDNA read counts. There are plenty of existing attempts to solve the issues by adjusting the assumed models, including research by the authors of commonly-used bulk differential expression tools (e.g. here).


You can use gffread to convert gff to gtf2, below is from the manual:

  In order to see the GTF2 version of the same transcripts, the -T
  option should be added:

gffread -E annotation.gff -T -o- | more

  The examples above also show that gffread can be used to convert a
  file between GTF2 and GFF3 file formats.


Changing the name of elements in a column

If your data is in a dataframe you can make use of plyr. Using mapvalues() from this library, you can rename provide a list of keys and a list of values. Instead of explicitly listing the keys, we can just use unique(mat$Tumor_Sample_Barcode) [credit to @dc37]. I use the paste0 function here and the colon operator to generate the vector of A1,A2,...,A21.

Assign the vector produced by mapvalues() to the column and you're off!

install.packages("plyr")
mat$Tumor_Sample_Barcode <-
plyr::mapvalues(x=mat$Tumor_Sample_Barcode,
                from=unique(mat$Tumor_Sample_Barcode),
                to=paste0("A", 1:21))

I'd recommend keeping this column and adding a new column with your simplified sample names, but this is not what you requested.


Given a Genomic Ranges of SNPs, how to inject these SNPs in genome via BSGenome?

This seems relatively complicated given the structure of a BSGenome object.

The creator of the package answered this question previously on the Bioconductor support forums:
https://support.bioconductor.org/p/86665/#86757

  We don't provide an easy way to inject arbitrary SNPs in an arbitrary
  BSgenome at the moment. However, it should not be too hard to forge
  the BSgenome ... package. First you would need to
  compute the sequences of the mutated chromosomes (you can
  use replaceLetterAt for this), then write them to a 2bit file (put
  them in a DNAStringSet object and call rtracklayer::export on it),
  then use that 2bit file to forge the BSgenome ....
  package (see the BSgenomeForge vignette in the BSgenome package for
  how to do this).
  - Herve Pages

Here is another similar question with more info about replaceLetterAt: https://support.bioconductor.org/p/26199/