Is there a command line tool to split a SAM/BAM file by CB (cell barcode) tag? I have a BAM file from a single cell sequencing experiment. Each read has had the cell barcode annotated in the CB tag. Some reads do not have a CB tag. For each of an arbitrary set of barcodes, I would like to split out the reads for each of those barcodes into a separate file. Is there a tool that does this?

I wrote a python script based on pysam that is free for anyone to use:

"""
MIT License
Copyright (c) 2020 Warren W. Kretzschmar
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

from pathlib import Path

import pysam

class BamWriter:
    def __init__(self, alignment, barcodes, prefix):
        self.alignment = alignment
        self.prefix = prefix
        self.barcodes = set(barcodes)
        self._out_files = {}

    def write_record_to_barcode(self, rec, barcode):
        if barcode not in self.barcodes:
            return
        if barcode not in self._out_files:
            self._open_file_for_barcode(barcode)
        self._out_files[barcode].write(rec)

    def _open_file_for_barcode(self, barcode):
        self._out_files[barcode] = pysam.AlignmentFile(
            f"{self.prefix}_{barcode}.bam", "wb", template=self.alignment
        )

def main(input_bam, barcodes_file, output_prefix, contigs):
    """Split a 10x barcoded sequencing file into barcode-specific BAMs

    input:
    barcodes_file: a file containing barcodes or a single barcode
    contigs: '.' for all contigs, 'chr1' for the contig 'chr1',
    or '1-5' for chromosomes 1, 2, 3, 4, and 5
    """
    alignment = pysam.AlignmentFile(input_bam)
    if Path(barcodes_file).is_file():
        with open(barcodes_file, "r") as fh:
            barcodes = [l.rstrip() for l in fh.readlines()]
    else:
        barcodes = [barcodes_file]
        print(f"Extracting single barcode: {barcodes}")
    writer = BamWriter(alignment=alignment, barcodes=barcodes, prefix=output_prefix)
    if contigs == ".":
        print("Extracting reads from all contigs")
        recs = alignment.fetch()
    else:
        if "-" in contigs:
            start, end = contigs.split("-")
            print(f"Extracting reads from contigs {start} to {end}")
            recs = (alignment.fetch(str(contig)) for contig in range(start, end + 1))
        elif "," in contigs:
            contigs = contigs.split(",")
            print(f"Extracting reads from contigs {contigs}")
            recs = (alignment.fetch(str(contig)) for contig in contigs)
        else:
            print("Extracting reads for one contig: {contigs}")
            recs = (alignment.fetch(c) for c in [contigs])

    for rec in recs:
        try:
            barcode = rec.get_tag("CB")
            writer.write_record_to_barcode(rec=rec, barcode=barcode)
        except KeyError:
            pass

if __name__ == "__main__":
    import sys

    main(
        input_bam=sys.argv[1],
        barcodes_file=sys.argv[2],
        output_prefix=sys.argv[3],
        contigs=sys.argv[4],
    )
```


What is the difference between NGS, 2GS, SBS and HTS? I've come across a bit of confusion about the initialism NGS, so think it would be a good idea to clarify this term (and similar terms like 2GS, SBS, and HTS) for this site with a bit of discussion.

What are the most common initialisms used to describe different sequencing technologies?

Here are my attempts at definitions:

Sanger: A method of sequencing that depends on chain-terminatiing dideoxynucleotides. This sequencing uses the differential flow of DNA sequences of different lengths through a gel to determine the original DNA sequence, producing a single sequence per reaction container.

NGS: Next-generation sequencing, also referred to as 2GS (second-generation sequencing). This term is used to describe the first wave of sequencing technologies that followed Sanger sequencing technology. The use of NGS has become more confusing with the advent of long-read sequencing, because it's a common assumption that "next-generation" refers to the most recent technology (which is incorrect in this case).

HTS: High-throughput sequencing. This term describes any type of sequencing technology that produces large amounts of data, usually in the form of millions of different sequences produced from the same sequencing run.

SBS: Sequencing by synthesis. This term describes a method of sequencing that depends on the synthesis of [DNA] bases in order for sequencing to be carried out. This definition can extend into long-read sequencing (e.g. PacBio sequencers depend on synthesis during sequencing), but is more typically associated with only the second-generation sequencing technology.


While opening a Bam file with the SeqAn library I get 'seqan::FileOpenError' When I try to open a simple bam file with the SeqAn library using the following code:

CharString bamFileInName = getAbsolutePath("/home/omair/TFG/Files/9827_2#49.bam")
BamFileIn bamFileIn(toCString(bamFileInName));

I get the following error:

terminate called after throwing an instance of 'seqan::FileOpenError'
  what():  Could not open file /home/omair/Programs/seqan//home/omair/TFG/Files/9827_2#49.bam: iostream error

Process finished with exit code 134 (interrupted by signal 6: SIGABRT)

The file exists and the path is correct.

Remove the getAbsolutePath call — your path is already absolute. getAbsolutePath garbles it, as you can see in the error message.


API/website for blasting peptide sequence against database of all plant proteomes I would like to blast a peptide sequence against all available plant proteome databases, i.e. a blastp metaproteome analysis.

What tools are available? Which database can I use?

A simple BLAST on the UniProt website does the trick (target database plants).

Best matches are subsequences ETVEREKE and EKKRRETV in Green alga and Eelgrass respectively.


Visualize PLINK pedigree files Is there a pedigree viewer that works with PLINK pedigree files? I saw an old post by Pierre, but I don't see the script on his blog.

We use kinship2 R package to visualize pedigree files (.fam plink files). For some reasons, some characters from plink ped aren't supported (like "-9" for unknown genre or affected status).
Here is the R script we use:

args<-commandArgs(TRUE)
pedPath<-args[1]
pdfFile<-args[2]

library(kinship2)
pedfile <- read.table(pedPath,sep="\t")
pedfile$V5 <- apply(pedfile,1,function(data){if(data[5] == 1 || data[5] ==2){data[5]}else{3}})
pedfile$affected <- apply(pedfile,1,function(data){if(as.numeric(data[6]) == 1){res <- data[6]}else{if(as.numeric(data[6]) == 2){res <- data[6]}else{res <- NA}}})

familyIDs <- as.character(unique(pedfile$V1))

pdf(file=pdfFile)
for(famID in familyIDs){
    ped <- pedigree(id=as.character(pedfile$V2[pedfile$V1==famID]),dadid=as.character(pedfile$V3[pedfile$V1==famID]),momid=as.character(pedfile$V4[pedfile$V1==famID]),sex=as.numeric(pedfile$V5[pedfile$V1==famID]),affected=as.numeric(pedfile$affected[pedfile$V1==famID]))
    plot.pedigree(ped,cex=0.7,align=TRUE)
    title(famID)
}

And we call it this way:

Rscript <kinshipScript.R> <plink ped file> <pdf output>


Are there any databases for gene co-expression or expression pattern clustering? I am currently working on gene clustering based on co-expression pattern in mouse brain. The problem is I do not have some solid way to test my result. Are there any suggestions for databases indicating groups of genes that co-express (in mouse brain), or even better indicating some similarity metrics between each pair of genes based on expression pattern?
Also, I do want some databases for differentially-expressed genes ;)
PS: By "differentially-expressed" I am referring to a pair of genes whose expression patterns in a single individual differ.

How about:

http://www.informatics.jax.org/expression.shtml
http://www.emouseatlas.org/emage/home.php

Other potential databases of interest:

https://biokeanos.com/search?q=co-expression+%2Bmouse


mrbayes, Sumt ; Error when setting parameter "Filename" (1) question 1
I wanted to use sumt command after completion of mcmc run. Whenever I use, sumt Filename=infile.nex, it gives me the error:
A matrix must be specified before sumt can be used
Error when setting parameter "Filename" (1)

But for the same analysis, I can use sump without any problem. What should I do now?
question 2
Another issue, I used burn in as 25% where as in Tracer, it is always showing 10%. No matter how much burn in I use, Tracer always shows 10% as burn in.
system stats
Operating system : Centos7
Version of MrBayes: 3.2.6

Thanks in advance,
Best Regards,
Zillur

Having thought this through, I think you've made a mistake. The file it wants is the outfile, you have given it the infile only. My original answer is below.

====================

From memory (my personal memory, not RAM, SSD or Stata), in the run command, you need to specify,

nruns=2

sumt will then combine both runs. If the file name is mystuff.run and nruns=2, there will be an output of mystuff.run1.t and mystuff.run2.t.

If nruns isn't set then sumt will not work, or if the MCMCMC output ends up with different names again sumt will not work.

With the burnin, you simply specify this in MrBayes. It is easy to write a script that will perform a 25% burnin. Tracer is nicely written, presumably C, its fast, but it has been around for a long time.


Modifying protein_scale in ProteinAnalysis in SeqUtil in Biopython? Is there a way to give weight arguments to each of the residues when using protein_scale in the ProteinAnalysis module of Bio.SeqUtils.ProtParam? I thought this would be done with the window and edge arguments, but from the documentation and the question How to allow "half windows" in biopython? I understand that the weights are symmetric and calculated automatically given a window size.
I'd like to understand what the idea behind this was, just to check that I understand this correctly. For example, if one calculates eg. hydrophobicity of a protein, why should it be a value per residue that depends on a particular window?
I was looking for a functionality like where the value given for each residue in the ParamData dictionary would be multiplied by a specified weight, like so:
```Y = ProteinAnalysis('MSTAEFAQLL')
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10)
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,edge[0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.])```

Of course I can code this myself, but was wondering if that functionality already existed.

There is not a way to control the edge weights as written, but making a modified version is straightforward. Namely make a class that inherits the class Bio.SeqUtils.ProtParam.ProteinAnalysis and overwrites the method relevant methods.
Having a gander at the files one sees the following:
def protein_scale(self, param_dict, window, edge_list: list):
    """[...]
    """
    # generate the weights
    #   _weight_list returns only one tail. If the list should be
    #   [0.4,0.7,1.0,0.7,0.4] what you actually get from _weights_list
    #   is [0.4,0.7]. The correct calculation is done in the loop.
    weights = self._weight_list(window, edge)
    scores = []
    ...

So the weights are generated by the 'private' method ._weight_list. Underscore means it is for internal use and for production code it is bad to assume that the name of the private method many not change, but Biopython has not changed in years so you'd be safe. And this is for a quick hack anyway...
class JankyProteinAnalysis(Bio.SeqUtils.ProtParam.ProteinAnalysis):
    overriding_weights = None

    def _weight_list(self, window, edge):
          if self.overriding_weights is not None: # dodgy override!
               print("Weight method hack done!")
               return overriding_weights
          else: # behave normally
               return super()._weight_list(window, edge)

So you'd do:
Y = JankyProteinAnalysis('MSTAEFAQLL')
Y.overriding_weights = [0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.]
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,float('nan'))

Note that the edge argument is no longer doing anything.


Is there any tools that can generating FBA model from genome sequences or proteome sequences? Suppose I have all sequences of genes or proteins of certain species, are there any free tools I can use to generate a flux balance analysis model out of them?

What I was thinking is more like to annotate all the enzymes based on some annotation method. Then mapping these enzymes to known meta-metabolic networks that we know for instance the meta-metabolic network in KEGG database.

To generate flux balance, one should now the stoichiometric numbers of all the reactions that are taking place, as well as the conditions (that could affect the reactions). I am not sure we have enough knowledge, time and computationally power to do this.

At the moment I'm only aware of a paper simulating the Mycoplasma genitalium cell. But I'm not sure it was worth the effort as I haven't seen replicated or redone with other species, perhaps the model was not accurate enough to be worth it.

The tools used are here. The licenses are not a standard, but at least they seem free.


How to visualize called narrowPeak files in UCSC Genome browser or IGV? I have called peaks using MACS2. Then I got a narrowPeak file like this.

1   4491713 4491913 testONE_peak_1  37  .   2.86003 3.72580 1.66717 80
1   4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
1   4571554 4572026 testONE_peak_3  72  .   4.16005 7.21000 4.85200 282
1   4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
1   4784848 4785061 testONE_peak_5  31  .   2.24218 3.19160 1.21569 133

Then I tried to visualize it in UCSC Genome browser and I got following error.

Error File 'testONE_peaks.narrowPeak' - Unrecognized format line 1 of file: 1 4491713 4491913 testONE_peak_1  37  .   2.86003 3.72580 1.66717 80 (note: chrom names are case sensitive, e.g.: correct: 'chr1', incorrect: 'Chr1', incorrect: '1')

So I added chr to the chromosome column instead of 1 (Chromosome number).

chr1    4491713 4491913 testONE_peak_1  37  .   2.86003 3.7258  1.66717 80
chr1    4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
chr1    4571554 4572026 testONE_peak_3  72  .   4.16005 7.21    4.852   282
chr1    4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
chr1    4784848 4785061 testONE_peak_5  31  .   2.24218 3.1916  1.21569 133

When I tried to load the file again I got following error.

Error File 'testONE_peaks_NEW2.narrowPeak' - Error line 1 of custom track: thickStart out of range (chromStart to chromEnd, or 0 if no CDS)

Then I added following line to the top and I could load the file to UCSC Genome browser.

track type=narrowPeak name="Somite narrowPeak"
chr1    4491713 4491913 testONE_peak_1  37  .   2.86003 3.7258  1.66717 80
chr1    4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
chr1    4571554 4572026 testONE_peak_3  72  .   4.16005 7.21    4.852   282
chr1    4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
chr1    4784848 4785061 testONE_peak_5  31  .   2.24218 3.1916  1.21569 133

But the problem is I can not see anything in the Genome browser. Only the name is there. But nothing else. Do you know what I should do?

I think you might have changed the separator (or at least have some kind of inconsistency from the required format) for your file. Note that peak output files from MACS2 are variants of BED files.

It seems you need to have a tab separator for this type of file.

I copied your example file and then ran awk -v OFS="\t" '$1=$1' your_peaks.narrowPeaks > fixed_your_peaks.narrowPeaks on it. This produced a result that I was able to see in both the UCSC genome browser and IGV.


Single-sample vs. joint genotyping I am trying to understand the benefits of joint genotyping and would be grateful if someone could provide an argument (ideally mathematically) that would clearly demonstrate the benefit of joint vs. single-sample genotyping.

This is what I've gathered from other resources (Biostars, GATK forums, etc.)

Joint-genotyping helps control FDR because errors from individually genotyped samples are added up, and amplified when merging call-sets (by Heng Li on https://www.biostars.org/p/10926/)

If someone understands this, can you please clarify what is the difference on the overall FDR rate between the two scenarios (again, with an example ideally)

Greater sensitivity for low-frequency variants - By sharing information across all samples, joint calling makes it possible to “rescue” genotype calls at sites where a carrier has low coverage but other samples within the call set have a confident variant at that location. (from https://software.broadinstitute.org/gatk/documentation/article.php?id=4150)

I don't understand how the presence of a confidently called variant at the same locus in another individual can affect the genotyping of an individual with low coverage. Is there some valid argument that allows one to consider reads from another person as evidence of a particular variant in a third person? What are the assumptions for such an argument? What if that person is from a different population with entirely different allele frequencies for that variant?

Having read several of the papers (or method descriptions) that describe the latest haplotype-aware SNP calling methods (HaplotypeCaller, freebayes, Platypus) the overall framework seems to be:

Establish a prior on the allele frequency distribution at a site of interest using one (or combination) of: non-informative prior, population genetics model-based prior like Wright Fisher, prior based on established variation patterns like dbSNP, ExAC, or gnomAD.

Build a list of plausible haplotypes in a region around the locus of interest using local assembly.

Select haplotype with highest likelihood based on prior and reads data and infer the locus genotype accordingly.

At which point(s) in the above procedure can information between samples be shared or pooled? Should one not trust the AFS from a large-scale resource like gnomAD much more than the distribution obtained from other samples that are nominally party of the same "cohort" but may have little to do with each other because of different ancestry, for example?

I really want to understand the justifications and benefits offered by multi-sample genotyping and would appreciate your insights.

Say you are sequencing to 2X coverage. Suppose at a site, sample S has one reference base and one alternate base. It is hard to tell if this is a sequencing error or a heterozygote. Now suppose you have 1000 other samples, all at 2X read depth. One of them has two ALT bases; 10 of them have one REF and one ALT. It is usually improbable that all these samples have the same sequencing error. Then you can assert sample S has a het. Multi-sample calling helps to increase the sensitivity of not so rare SNPs. Note that what matters here is the assumption of error independency. Ancestry only has a tiny indirect effect.

Multi-sample calling penalizes very rare SNPs, in particular singletons. When you care about variants only, this is for good. Naively combining single-sample calls yields a higher error rate. Multi-sample calling also helps variant filtering at a later stage. For example, for a sample sequenced to 30X coverage, you would not know if a site at 45X depth is caused by a potential CNV/mismapping or by statistical fluctuation. When you see 1000 30X samples at 45X depth, you can easily know you are looking at a CNV/systematic mismapping. Multiple samples enhance most statistical signals.

Older methods pool all BAMs when calling variants. This is necessary because a single low-coverage sample does not have enough data to recover hidden INDELs. However, this strategy is not that easy to massively parallelized; adding a new sample triggers re-calling, which is very expensive as well. As we are mostly doing high-coverage sequencing these days, the old problem with INDEL calling does not matter now. GATK has this new single-sample calling pipeline where you combine per-sample gVCFs at a later stage. Such sample combining strategy is perhaps the only sensible solution when you are dealing with 100k samples.

The so-called haplotype based variant calling is a separate question. This type of approach helps to call INDELs, but is not of much relevance to multi-sample calling. Also, of the three variant callers in your question, only GATK (and Scalpel which you have not mentioned) use assembly at large. Freebayes does not. Platypus does but only to a limited extent and does not work well in practice.

I guess what you really want to talk about is imputation based calling. This approach further improves sensitivity with LD. With enough samples, you can measure the LD between two positions. Suppose at position 1000, you see one REF read and no ALT reads; at position 1500, you see one REF read and two ALT reads. You would not call any SNPs at position 1000 even given multiple samples. However, when you know the two positions are strongly linked and the dominant haplotypes are REF-REF and ALT-ALT, you know the sample under investigation is likely to have a missing ALT allele. LD transfers signals across sites and enhances the power to make correct genotyping calls. Nonetheless, as we are mostly doing high-coverage sequencing nowadays, imputation based methods only have a minor effect and are rarely applied.


How do I generate a color-coded tanglegram? I want to compare two phylogenies and colour the association lines based on some metadata I have. I have been using ape cophyloplot but I have not had any success in getting the lines to colour accurately according to my data (see previous question).

Note that in my actual work flow I define the colour scheme using a palette to control the colour outcome.

I want a means to make a tangle using phylogenies which I can format. Preferably in R. I like to get an output like this:

I think you can try dendextend, in this manual there is an example of coloring the branches. I don't think it is exactly like your coloring, but with a little tweaking you might get your colorscheme in there.

The manual mentions an argument called color_lines for the function  tanglegram():

# The `which` parameter allows us to pick the elements in the list to compare
iris_dendlist %>% dendlist(which = c(1,4)) %>% ladderize %>%
   # untangle(method = "step1side", k_seq = 3:20) %>%
   set("rank_branches") %>%
   tanglegram(common_subtrees_color_branches = TRUE)

According to the manual, the code above produces an image like this:

Try to see if you can fit your color palette in there.


How to automate NCBI genome download I need to download all the completely assembled cyanobacterial genome's GenBank file(.gbff) from NCBI(RefSeq or INSDC ftp data).

For this I think, the steps are:

Need to find the completely assembled genomes.
find the GenBank file URL based on the taxonomic name.
download the GenBank file(.gbff file)

Is their any way to do this by using any python module or any other idea??

Outline of solution:

get this file: ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt
filter for lines having "Complete Genome" in column 12
filter for lines having a taxid (column 6) that is a descendant of taxon id 1117 (phylum Cyanobacteria), you can use nodes.dmp from ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
The file download URL is the concatenation of column 20 and the last field of column 20, after separating it at '/', followed by appending '_genomic.gbff.gz'.  For example: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/007/365/GCF_000007365.1_ASM736v1/GCF_000007365.1_ASM736v1_genomic.gbff.gz

Edit:
Since I will need something like that soon, I made a Perl script for downloading genomes.
For Cyanobacteria, you would do ./download-refseq-genomes.pl 1117


creating a tab delimited file I am working on a project using a fasta file. I am writing my command in nano within command-line and executing using python, also within my command-line.
I would like my command to provide me with a tab delimited file with three columns: first column should contain my sequence name, second column should provide me with my sequence length, and the third column should show the sequence itself.
I have written the following command so far within nano:
from Bio import SeqIO
import sys
for hello_fasta in SeqIO.parse(sys.argv[1], "fasta"):

  list = hello_fasta.split("\t")

  print hello_fasta.description
  print (len(hello_fasta.seq))

For example, I would like my command to provide me with the desired output and with the following order: Gene name ; Gene length ; Gene seq
H0192X 26 FORUWOHRPPTRWFAWWEAKJNFWEJ

You can use a list and insert() to add an element in a specific order, then expand the list with *. Or you can use join().
from Bio import SeqIO
import sys

for hello_fasta in SeqIO.parse(sys.argv[1], "fasta"):
  sequences = []
  sequences.insert(0, hello_fasta.description)
  sequences.insert(1, len(hello_fasta.seq))
  sequences.insert(2, hello_fasta.seq)
  # option 1
  print(*sequences, sep='\t')
  # option 2
  print('\t'.join(map(str, sequences)))


Where can I find summary data for how common certain mutation *types* are? I'd like to know how common certain mutation types are in public data sets like the 1000 Genomes, ExAC, and ESP6500. Specifically, I'd like to know the distribution of stop-gains, stop-losses, frameshift, and other mutation types. For example, what is the median number of stop-gains observed across individuals in the ExAC population? I'd like to have the full distribution, but median values would work.

To be clear, I am not looking for allele frequencies.

I didn't find this kind of summary data when I looked at the respective websites. Is there a good source for these summary data?

This dbsnp github repo seems to offer that information.

https://github.com/ncbi/dbsnp/tree/master/Build%20Announcements/151


are GSEA and other geneset enrichment analysis supposed to yield extremely different results between them? I have recently ran in R four geneset enrichment analysis in the same database (TCGA:breast cancer) comparing two intrinsic subtypes. The methods I used were:

  MIGSA, that imports mGSZ package and combines it with a SEA algorithm. Using RNA-seq TMM normalized counts as input, because the software supports it. All the parameters are as default.
  mGSZ with standard parameters, using RNA-seq TMM normalized log2 cpm. All parameters as default
  FGSEA preranked, using log fold change as ranking metric, obtained from limma deg from RNA-seq TMM normalized counts. All the parameters are as default.
  GSVA + limma: extracted by-sample enrichment score from RNA-seq TMM normalized log2 cpm, then ran a differential expression analysis with geneset enrichment scores as if they were simple gene expression values. I know this is not a formal approach but I wanted to compare. All the parameters are as default.

The genesets tested are the full combined list of KEGG, GO:BP, GO:MF and GO:CC. In all cases I trimmed the genesets to genesets of length <500, which leaves a total of 21723 genesets tested.
Sadly, results are confusing. For starters mGSZ for some reason is cropping the output and giving results for only around 9700 of those genesets -I've already sent an email to the package mainteiner asking about that. But that isn't the biggest problem, the most shocking thing is that none of the analyses performedt has a good correlation to the results of any of the other analyses. Here I plotted the -log10(adjusted p-value + 0.0001) of each analysis agaisnt every other analysis, with lines indicating -log10(0.05)=1.3 and -log10(0.01)=2.

I was expecting the results to be somehow different, but here I can't see even something remotely close to a correlation between the different results. It is also evident that MIGSA and GSVA+limma output a lot more significatively enriched genesets than the other two methods.
I have gone through the code time and again, and I can't find anything errors. Also I have searched for studies doing this kind of direct comparison between methods, without luck so far. Has anybody experienced this kind of inconsistencies between gene set enrichment analysis results with different tools? What could I possible be doing wrong?.

Each of these methods do something different, so it is reasonable to expect different results. The bottom line is that there isn't a single question you can ask when you do an "enrichment test". For example: you can test if a set of genes are more expressed than this other list of genes or they are above the mean levels of any other gene.
Even more, you use the gene ontology, which has an underlying structure (directed acyclic graph) that no method you used takes into account but could modify the results if taken into account. However, the KEGG pathways do not have this structure. So when you want to apply a method mixing both data sources you need to decide if you want to use it or not.
This is also complicated because for most of these methods you can only compare the p-value which depends on the underlying assumptions of the methods (uniform, normal, binomial distribution?)
The latest comprehensive study comparing enrichment methods including some of the methods you used and more is this one (that I know of): Toward a gold standard for benchmarking gene set enrichment analysis.
Last, what would a correlation tell you here? It is know that there is little agreement between methods. You couldn't say that one method is equal to the other unless you tested with several datasets and you had always a correspondence on the p-value.


how to get FlyBase ID conversion table FlyBase uses systematic naming of genes, transcripts and proteins in a format FBtypeNUMBER. However, the numbers are not the same for genes and their products. For example gene FBgn0259101 transcribes to FBtr0299513 and translates to FBpp0288787.

On this webpage users can conver individual IDs, but what If I need to convert practically whole transcriptome? Where can I get the full translation table?

Indeed the web interface does not seem to link the full translation table, but it is available on their ftp server at releases -> release_ID -> precomputed_files -> genes -> fbgn_fbtr_fbpp_fb_2019_06.tsv.gz

For example, this is the latest (FB2019_06) translation table available.


Running snakemake with multiple conda environments without rebuilding environments every time I have some snakemake pipelines that have rules which run in different environments. Is it possible to run the pipeline without having to rebuild the environments from scratch every time when run in a different directory?

This is the way I would normally run snakemake:

snakemake --directory <my_output_dir> --use-conda

My snakefile might have rules that need to use different environments:

rule A:
   ....
   conda: envA.yml
   ....

rule B:
   ....
   conda: envB.yml
   ....

So if I wanted to do two different analyses with the same pipeline, is there anyway that I can point them to the same environments without having to install everything separately in each output directory?

The snakemake documentation makes it seem like they don't currently support this. You can either use 1) --conda-prefix when running snakemake and then point to a single env or 2) use the --use-conda directive and the env.yml files by rule, in which case it will create new envs any time you run it in a new output directory. So my case doesn't seem to be covered, but does anyone have a workaround to suggest? I often find conflicts between some of the R packages I use and python packages, so I tend to keep them in their own environments to avoid headache.

I'm starting to get neck-deep in conda environments at this point which is taking up space from my disk quota, but also I am bothered by the fact that I may actually want to run analysis of data from different projects with the same environments. If the environments are rebuilt at each execution, I do not know that they are the same, short of specifying the exact version info for everything I use in the env.yaml files (and even then, they could be slightly different if the repo has changed).

I have looked at this similar question: Running Snakemake in one single conda env
in which it is suggested there is a way to do such a thing, but I'm afraid I don't understand how to implement it without an example.

Thanks

Each YAML file will always create its own environment. The packages in each aren't downloaded and installed each time, if they're already installed elsewhere they're just hard-linked in. You'll want to use both --use-conda and --conda-prefix. The latter will point to the base directory in which all conda environments will be stored. This has the benefit that running the same workflow in different places won't result in a new env being created (by default snakemake will do so inside a .snakemake directory otherwise).


How to generate two fasta files from samtools phase output? samtools phase is an easy-to-use tool to phase variants from a bam file and a reference.

For a sequence that has two haplotypes, like this:

   position: 123456789
haplotype 0: GACCATATT
       SNPs:  **    *
haplotype 1: GTACATACT

The output looks like this:

CC  M?   chr    PS     pos  allele0 allele1 hetIndex #supports0 #errors0 #supp1 #err1
CC
CC
PS      frame-1   2    9
M1      frame-1  60    2      A       T       1       3       0       23      0
M1      frame-1  60    3      C       A       2       35      1       24      0
M1      frame-1  60    8      T       C       3       46      0       30      0

How can I take the reference fasta file and the output of samtools phase to generate one fasta file for each of the haplotypes?

Answer January 20th, 2020. If you run samtools phase with the option -b my_prefix it outputs a bam file for both haplotypes, named my_prefix.0.bam and my_prefix.1.bam. If the sequence has been correctly phased, use these bam files to correct the reference using pilon. This is preferable to the method described in the January 14th, 2020 answer because pilon will correct SNPs and indels, whereas the python program below will only correct the specific SNPs that samtools phase found. samtools phase does not find indels.

Answer January 14th, 2020. This python program below takes a fasta file as an argument and the output of samtools phase piped in, then prints out a fasta file of the two haplotypes.

Run the program like this:

samtools phase reads.bam | phase2fasta.py ref.fa > haplotypes.fasta.

Here is the code of phase2fasta.py:

#!/usr/bin/env python

import os
import sys

if (len(sys.argv) == 1) or not os.path.exists(sys.argv[1]):
    raise Exception("The first and only argument is the reference fasta.")
allele0 = []
seqname = ""
with open(sys.argv[1]) as f:
    for line in f:
        if line[0] == ">":
            seqname = line[1:].split(" ")[0].strip()
        else:
            for c in line.strip().upper():
                allele0.append(c)
allele1 = [x for x in allele0]
for line in sys.stdin:
    if line[0:2] == "M1":
        fields = line.strip().split()
        if fields[1] == seqname:
            pos = int(fields[3]) - 1
            allele0[pos] = fields[4]
            allele1[pos] = fields[5]
print(">{}_allele0".format(seqname))
print("".join(allele0))
print(">{}_allele1".format(seqname))
print("".join(allele1))


Group genes by functional categories suming expression values Using the count of rpkm values from genes in a metagenome sample, I want to group these genes into established categories (for example KEGG or COG). For each sample, my goal is to determine which categories are better represented in each one.

Considering the example table whit genes, samples (1-5) and categories.
Each gene has their expression and the column categories classify them.

import sys
import numpy

dictionary = {}

table = open(sys.argv[1], "r")

for line in table:
    next()
    cols = line.strip().split("\t")
    gene = cols[0]
    s1 = float(cols[1])
    s2 = float(cols[2])
    s3 = float(cols[3])
    s4 = float(cols[4])
    s5 = float(cols[5])
    counts = [s1, s2, s3, s4, s5]
    kegg = cols[6].strip().split(" ")[0].replace('"','')

    if kegg not in dictionary:
        dictionary[kegg] = [counts]
    else:
        dictionary[kegg].append(counts)

#print dictionary
for k, v in sorted(dictionary.iteritems()):
    m = numpy.array(v)
    print k,  [sum(m[:,i]) for i in range(5)] # range(number_of_samples)

Original table to count the columns values:

gene    s1  s2  s3  s4  s5  Category
name01  0   2   2   0   0   A
name02  3   0   1   0   0   A
name03  0   2   1   0   0   A
name04  0   0   1   0   0   B
name05  5   0   1   0   0   C
name06  1   0   0   0   0   D
name07  2   0   0   0   0   D
name08  0   0   3   0   0   E
name09  1   0   0   0   0   F
name10  1   0   0   0   0   F
name11  3   0   0   0   0   F

The result obtained is this:

"""
# eliminated from the question

Cat s1  s2  s3  s4  s5
A   0.0 2.0 1.0 0.0 0.0
B   0.0 0.0 1.0 0.0 0.0
C   5.0 0.0 1.0 0.0 0.0
D   2.0 0.0 0.0 0.0 0.0
E   0.0 0.0 3.0 0.0 0.0
F   3.0 0.0 0.0 0.0 0.0
"""

But my expected (obtained) table is this:

Cat s1  s2  s3  s4  s5
A   3.0 4.0 4.0 0.0 0.0
B   0.0 0.0 1.0 0.0 0.0
C   5.0 0.0 1.0 0.0 0.0
D   3.0 0.0 0.0 0.0 0.0
E   0.0 0.0 3.0 0.0 0.0
F   5.0 0.0 0.0 0.0 0.0

As suggested in the comments, the pandas module is quite convenient for this type of work.

Here is how you could do (assuming you have tab-separated values as input, and you want similar file format for the output):

import pandas as pd

counts = pd.read_table("counts.tsv")
# At this point, counts looks as follows:
#       gene  s1  s2  s3  s4  s5 Category
# 0   name01   0   2   2   0   0        A
# 1   name02   3   0   1   0   0        A
# 2   name03   0   2   1   0   0        A
# 3   name04   0   0   1   0   0        B
# 4   name05   5   0   1   0   0        C
# 5   name06   1   0   0   0   0        D
# 6   name07   2   0   0   0   0        D
# 7   name08   0   0   3   0   0        E
# 8   name09   1   0   0   0   0        F
# 9   name10   1   0   0   0   0        F
# 10  name11   3   0   0   0   0        F

summed_by_cat = counts.groupby("Category").sum()
# At this point, summed_by_cats looks as follows:
#           s1  s2  s3  s4  s5
# Category
# A          3   4   4   0   0
# B          0   0   1   0   0
# C          5   0   1   0   0
# D          3   0   0   0   0
# E          0   0   3   0   0
# F          5   0   0   0   0

# Change the index name to get the desired first column header in the output:
summed_by_cat.index.name = "Cat"

# Write to a file
summed_by_cat.to_csv("by_cat.tsv", sep="\t")

And resulting the file contains the following:

Cat s1  s2  s3  s4  s5
A   3   4   4   0   0
B   0   0   1   0   0
C   5   0   1   0   0
D   3   0   0   0   0
E   0   0   3   0   0
F   5   0   0   0   0

You could also load your table using counts = pd.read_table("data.tsv", index_col="gene") and counts would have the gene names as index instead of default numbers:

        s1  s2  s3  s4  s5 Category
gene
name01   0   2   2   0   0        A
name02   3   0   1   0   0        A
name03   0   2   1   0   0        A
name04   0   0   1   0   0        B
name05   5   0   1   0   0        C
name06   1   0   0   0   0        D
name07   2   0   0   0   0        D
name08   0   0   3   0   0        E
name09   1   0   0   0   0        F
name10   1   0   0   0   0        F
name11   3   0   0   0   0        F

This doesn't seem to affect the final outcome.


How to compute Shanon entropy I want to analyze some viral sequences. I want to use shanon entropy as a measure to analyze the sequences. However I do  not understand the concept of shannon entropy. Would anyone explain please how it works or advise on some resources i can use to understand shannon entropy.

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

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

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

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


How to calculate p-values for fold changes? I'm currently working with data from a Luminex multiplex assay. In this assay, the concentrations of 17 different analyte proteins were identified in 12 groups in triplicate. One of these 17 groups was used as the control, and the log2 fold changes were calculated for the analyte concentration of each sample in each group using the average control concentration for that analyte.

However, now I would like to calculate a p-value for the identified fold changes if possible.

My current preliminary idea is to perform the t test for each group compared with the control group for each analyte, but I don't think that my group sizes are large enough for this. My statistical knowledge is lacking, so if there is a better way of calculating these p-values I don't know about it.

Is there a better way of obtaining these p-values?

Your null hypothesis would be that the fold-changes are 0, so you can either do the T-test accordingly or simply do away with the fold changes and perform the T-test between the raw values in the group (this is preferable to performing a T-test of one group vs. 0, since it allows you to assess the expected variability around 0). Note, however, that 3 samples is pretty much the bare minimum needed for any statistics, so your power will be terrible and your estimation of variance will likely be pretty inaccurate. Consequently, take any p-values with an appropriate grain of salt.

To my knowledge there aren't any good alternatives for your situation without having a good expected background distribution.


Cuffmerge: EOF marker is absent. The input is probably truncated I need to merge my all transcripts.gtf from cufflinks output follow this command line :
cuffmerge -o merged_gtf_output -p 15 -s ref.fasta -g anot.gtf assembly.txt

But I have this warning :
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File merge_B3//tmp/mergeSam_fileOaNffV doesn't appear to be a valid BAM file, trying SAM...
[21:40:04] Loading reference annotation.
[21:40:05] Inspecting reads and determining fragment length distribution.
Processed 3921 loci.
> Map Properties:
>       Normalized Map Mass: 33242.00
>       Raw Map Mass: 33242.00
>       Fragment Length Distribution: Truncated Gaussian (default)
>                     Default Mean: 200
>                  Default Std Dev: 80

Can someone tell me what this warning means and how to fix it?

There is a post about this on seqanswers. The poster says the following about cuffmerge behavior:

  I kept the temporary files while running cuffmerge. It builgs a SAM file from transcripts.gtf files. Then it seems first it reads it as BAM since I read samtools take the files as BAM default. Then it finds it is not a BAM so it tries SAM. Then it can read it and resume the task.

So it sounds like the BAM file is generated from GTF as a kludge inside cuffmerge, not the other way around.

Did you change any of the intermediate files or delete anything? Is there maybe a different cuffmerge run that it might be colliding with?

Did you get any output? it sounds like this error is frequently observed when actually the command completes fine, e.g. https://www.biostars.org/p/189139/. See other posts by googling the problem maybe, e.g. https://www.biostars.org/p/344439/

(I know this isn't a proper answer but it was too involved for a comment.)


PDB file downloading: pymol automation vs. manual I automated a PDB download using a Pymol script (below)
python
pdb_lists = ['3LZG', '6HJN'] # lots more pdbs
for x in pdb_lists:
    cmd.fetch(x)
    cmd.select(x)
    cmd.save(x + '.pdb', x)
    cmd.delete(x)
cmd.quit()
python end

When the script is run within Pymol it pulls down a collection of pdbs, however this flopped when I  plugged it through a Biopython PDB parser, but was okay when I manually downloaded the same pdb file from RCSB.org
A manual downloaded PDB was perfect, meaning the amino acids which cannot be observed by X-ray crystalography (e.g. they not sufficiently stationary) were recognised as being absent from the parser (resulting in a protein ~390 amino acids), but the Pymol scripted download contained all amino acid residues, i.e. even those which couldn't be observed in the X-ray (resulting in a protein of ~440 amino acid residues). Any ideas why the script resulted in a pdb that flopped and how to correct it?

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.


How do I tell featureCounts to ignore fragments that overlap with more than one feature? I am using featureCounts from the subread package to count the number of fragments per gene in an RNA-seq experiment. I feed it a BAM file with paired-end reads aligned to a reference genome. In the help page of the executable, which I have installed using conda install subread, there is a section about overlaps, including for example the -O option which will assign the reads/fragments to all overlapping features. What is the default behaviour? Assigning to reads/fragments to the feature with largest overlap? not assigning multi-overlapping reads/fragments?

I would appreciate your help on this.

OS: macOS X

featureCounts version: 1.6.3

If a fragment overlaps multiple genomic features then it's unassignable, which is the appropriate way to handle ambiguous mappings when the resulting counts are then given to tools like DESeq2.


How to estimate the phenotypic variation explained by top SNPs from a GWAS study? I have conducted a large-scale GWAS study and got a few significantly associated SNPs. I used GEMMA with -lmm 1 options to run the GWAS and obtain the beta and standard-error estimates. I want to estimate the percent phenotypic variation explained by each of the significant SNPs. I used the following procedure for estimating the variance explained in R:
fit <- lm (Phenotypic_value ~ SNP_data, data = a)
summary(fit)$adj.r.squared

Here, the datafile a contains three columns namely, sample_ID, Phenotypic_value for each sample, and the biallelic SNP_data. I got a value which is 0.43 meaning 43% phenotypic variation explained by the SNP.
Again, I used another formula which is: 2*f*(1-f)*b.alt^2. Here, f is the minor allele frequency and b.alt is the effect size i.e. beta estimate obtained from GEMMA. This gives me a value of 0.03 meaning 3% variation explained which seems reasonable to me.
My question is that which of the above method is correct? Why do I obtain such a huge difference between the two approaches? Is there any other way to estimate the percent variation explained?
Alternatively, from the GEMMA google group, I have got this formula pve <- var(x) * (beta^2 + se^2)/var(y). But I do not understand how can I obtain the value of var(x) and var(y).
I borrowed the first and the last formula from GEMMA google group discussion. The first one is obviously a normal linear model equation in statistics. The second formula is from link. Also, the data structure is three columns, 1 is with individual, 2 is the genetic information (individual carrying A or T) and 3 is the phenotype.
In the last formula, var(x) is the variance of the genotype vector as they said in GEMMA group and the beta is the effect size estimate and se is the standard error from GEMMA output.
It will be great to receive some feedback on this. Thank you.

Answer from @maximilian-press, converted from comment:
At a guess, the GEMMA linear mixed model is controlling for population structure and other SNPs, which your R example is not (it only fits 1 SNP?), but it's hard to be sure with this information. It would be helpful to see data samples. There is probably more context and usage information available in those resources beyond just the equations. It is especially difficult to evaluate with no terms defined as in your last example (defining the terms in the earlier ones is a good start).


How are centromeric/telomeric regions categorized in RepeatMasker? "Other" category? I haven't been able to find the following in the documentation at http://www.repeatmasker.org/

(1) Why are certain regions classified as "Other"? Are these regions impossible to concretely classify due to alignment issues?

(2) I believe that centromeric and telomeric regions would be classified as tandem repeats. Are telomeric/centromeric regions classified as "DNA" within the hg19 version of RepeatMasker?

Are these documented somewhere?

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

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


if I know the number of sequencing circles can I give this information to DESeq2? I am trying to understand library normalization in DESeq2. I would like to ask the following: I know that some samples have been run 15 cycles and some others 20, can I give this information to DESeq2, would it be useful? I mean, I know that DESEQ2 uses an algorithm to normalize library size, but since I have this piece of information, I wonder if it would be useful to provide it somehow to DESEQ2. If it is not, would it make sense to incorporate this information in the design matrix, something like: design(~cycle+condition)? Thank you a lot!

You can include that info in the design, though you'd have to decide whether to include it as a number of just a factor.  I've never heard of anyone including PCR cycles.  I've heard of people including RIN scores, though I've never done it.
If the PCA looks alright, I wouldn't worry about.  I've got tons of examples of the lab people saying "I did a few more cycles on these few samples" and then you can't tell any difference in the PCA.  I think needing to do more cycles likely has more to do with variance in the prep than a real biological difference.


Smith-Waterman Pairwise Local Alignment Algorithm I have written a R code for Smith-Waterman (SW) algorithm.

I wanted to check the results of my code with online examples. In all cases -except one- I found the same alignment. So there is this one example -let's say example A*- for which I found a quite different alignment result: in the example they give three equivally good alignments whereas I found only one that corresponds to no none of those three alignments (we have the same score, though).

I am now trying to understand whether example A is wrong or if my code is not working properly.

Therefore, I have two requests:

Does anyone have a dataset of the result of SW used on two sequences?
I could then compare these results with the output of my code.
Does anyone have whatever way to compute SW algorithm? Then, I could give that person the example A, which could then try this example A, and give me the output.

This would be very helpful if someone could help, I have already spent so much time on my code, without finding anything that could explain those different outputs.

Thank you!

Example A: Sequence A = "TGCTCTGT" and Sequence B = "TACTATCAT"

There are many implementations of the Smith-Waterman algorithm. However, they are not all guaranteed to give the same output because although it's an optimal local alignment algorithm, it still depends on the chosen scoring scheme (i.e. substitution matrix and gap penalties).

With this in mind, you could try EMBOSS Water online:
https://www.ebi.ac.uk/Tools/psa/emboss_water/

Or in Biopython, the Align.PairwiseAligner() (which uses SW by default when aligner.mode is set to local):
https://biopython.org/DIST/docs/api/Bio.Align.PairwiseAligner-class.html


use same output in two processes in nextflow dsl2 This is my workflow:
pre_align()
pre_align.out.single_fastqs.view()
get_fq_info(pre_align.out.single_fastqs)

align_bwa(get_fq_info.out.fq_info)

align_bowtie2(get_fq_info.out.fq_info)

where I want to use the same output from get_fq_info as input for either align_bwa and align_bowtie2`.
As it's shown here only align_bwa is executed. If I comment align_bwa THEN align_bowtie2 gets executed.
How can I fix this?
Thank you!

I have solved this problem in the past by duplicating the output channel of the desired process in order to feed the same output into multiple downstream processes.


Putting labels of different sizes on one PyMOl Object I'm new to PyMol (and StackExchange!) and working on my first project. I have the structure of a protein as an object, called PolyA-M, and the idea is that the residues of it are shown as spheres of differing sizes based on a calcualted conservation value. I want to label each sphere with it's amino acid, but to have the label size corresponding to sphere size. Here is a sample of the code that might be used to label one amino acid.

alter ( resid 138 ), resn = "Q135"
alter ( name CB and resid 138 ), vdw = vdw * 0.8 * 0.679934640522875
set_color col138, [ 1.0, 0.181, 0.181 ]
show spheres, name CB and resi 138
show sticks, ( name CA or name CB ) and resid 138
set label_size, 14, PolyA-M
label PolyA-M and name CB and resid 138, "Q"
color col138, name CB and resi 138

In this case, residue 138 is medium sized sphere with a size 14 label Q that fits well. However, if later on in the code another residue is labelled as so:

alter ( resid 198 ), resn = "A208"
alter ( name CB and resid 198 ), vdw = vdw * 0.8 * 0.359803921568627
set_color col198, [ 1.0, 0.601, 0.601 ]
show spheres, name CB and resi 198
show sticks, ( name CA or name CB ) and resid 198
set label_size,  8, PolyA-M
label PolyA-M and name CB and resid 198, "A"
color col198, name CB and resi 198

This will label residue 198 with a small A to fit the small sphere. However, this will also change the earlier label to size 8 as well, making it too small for the sphere.

Is there any way for me to prevent this, such that each residue keeps its own label size?

I know that one way is to create multiple identical objects and keep all labels of the same size restricted to the same object (For example, PolyA-M_10 contains all residues with label size 10, polyA-M_14 all residues with label size 14) but I was wondeirng if there is a more efficient way?

Thank you for any help! This is my first Stack Exchange question so any feedback would be appreciated :)

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

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


significant SNPs to annotated candidate genes I have performed GWAS and got some significant SNPs. How can I get candidates genes in the flanking regions of those SNPs and also their functional annotations (eg. KEGG).

You can use BEDOPS' closest-features tool. You only need to have the SNPs coordinates in .bed format and a .bed file with gene locations. BEDOPS will look for the closest genes upstream and downstream of your SNPs. The command will look something like this:

closest-features yourSNPs.bed yourAnnotation.bed

Then, you can obtain your functional annotations through different tools such as biomart.


e-utils api: get chromosome position for each snp on a gene So I have the following e-utils api-link:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi?dbfrom=gene&db=snp&id=5726&retmode=json

it returns all the SNP-ids for the gene with the id 5726 as an example. Additionally to the SNP-ids I also need the SNP position and the genomic distance. Also I would need the chromosome but I could query that from the gene as well rather than for each SNP. Do you know how I can add all these informations to the resulting json/xml?
The first one in the resulting list is the SNP with the id 1490381980. Further information about it can be found on dbSNP: https://www.ncbi.nlm.nih.gov/snp/?term=1490381980%5BUID%5D
I really want to avoid creating a request for each single SNP and would prefer to query everything at once. Due to the unavoidable performance-issues. The Tas2R38 doesn't have that many SNPs but there are other genes with way more.

If your dataset is limited to human SNPs, and you know the list of SNPs, I would advise taking a look at the dedicated NCBI Variation Services.  This would be one request for each SNP, but if you fire off 1-2 per second, I think you'll be ok.
For example, for rs1490381980, you could try something like:
curl -X GET 'https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/1490381980' \
   | jq '.primary_snapshot_data.placements_with_allele[] | "\(.seq_id) \(.alleles[0].allele.spdi.position)"'

It would be better to query for the SNP and parse the JSON in python, but jq helps provide an idea of what is possible.
Alternatively, you could download the JSON or VCF data from the NCBI dbSNP FTP site, and process out the information you need.


HVG and variation coefficient I'm looking for highly variable genes from a data set. I have calculated the coefficient of variation of the expressed genes of my dataset and inserted them in a histogram graph. Now I should select only the highly variable genes. Should I select genes with a coefficient of variation >1, or how?
thank you for your help,

regards

There seems to be a dividing line around CoV=0.9.  Perhaps this is a good place to start for a cutoff.


Parallelizing WGCNA k-means clustering and merging smaller clusters WGCNA carefully takes advantage of server parallelization in several functions, and I can switch this on with:

enableWGCNAThreads(nThreads=28)

That speeds up some functions significantly, but the two bottleneck functions in the pipeline (in my experience) are part of the "Projective k-means" operation, primarily "..k-means clustering.." and "..merging smaller clusters..".

Unfortunately, parallelization does not work for those two functions.

Am I missing something?  How can I parallelize these functions?  Otherwise I'm wasting HPC resources and cluster time.

For reference, I'm working on a 28-processor HPC node with 96GB RAM, and I'm trying to run analysis on a 15,000-square correlation matrix.

The bad news is that, indeed, projectiveKMeans is not parallelized and I am not sure how much of it is (easily) parallelizable. The good news is that with 15k features (genes) and 96GB RAM you don't need to run preclustering at all. Just analyze the whole data set in one block.

If you for some reason insist on splitting the data into blocks, make sure the BLAS R uses on the cluster is a fast BLAS, ideally one that can parallelize matrix operations. A lot of time in merging of small clusters is spent calculating singular value decompositions, which becomes much faster using a good BLAS. HPC clusters will often have Intel KML available, or the admins could install OpenBLAS.


Optimal design for an RNA-seq experiment with ERCC RNA Spike-In Control Mixes I am planning an experiment on gene expression. I have 6 conditions and 3 replicates in each condition. I would like to use ERCC spike-in mixes for quality check: lower limit of detection, dynamic range and also check that the fold-changes are accurate. For example, by knowing that spike-in sequence xyz in condition 2 should be a 4-fold higher than in condition 1, does it actually get detected as a 4-fold change in my analysis? To assess the fold-changes, I need to use 2 different mixes. Only one mix should be used per sample. The manual suggests either adding mix 1 in a reference and mix 2 in all other samples or randomly assigning mixes 1 and 2 among samples. The latter is not very clear to me. Should I use the same mix for all replicates within a condition or should I randomise really among all the samples? Or can I randomise among conditions (for example 1, 3, 4 get mix 1 and 2, 5, 6 get mix 2)? So far, I am intending to do the latter, so using mix 1 in conditions 1, 3, 4 and mix 2 in conditions 2, 5, 6. This would also be cheaper because I can only buy equal volumes of mix 1 and 2. Do people have any advice on this?

I have 6 conditions and 3 replicates in each condition.

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

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

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

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

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

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


Is there a more elegant solution to the bwa-mem: paired reads have different names error? I'm currently trying to run bwa-mem on Influenza substrains using the following command:

~/bwa mem h5n1_1_cons.fa h5n1_1_read1.fq h5n1_1_read2.fq

h5n1_1_cons.fa is the consensus sequence for substrain h5n1_1, and the fq files are paired end reads 1 and 2. Running the program gives me this error:

[mem_sam_pe] paired reads have different names: "JQ714243-920-1", "JQ714243-920-2"

Visual inspection shows that the read names are identical until the last character. I have found a package (bbmap) with a utility, bbrename.sh, that renames the sequences, producing:

./bbrename.sh in=h5n1_1_read1.fq out=fixed prefix=h5n1_1

which yields the following output

-bash-4.2$ head fixh5n1_1_read1.fq
@h5n1_0
gaaaatgtctcaggagttgatgcatggtgtttaatctttggtttgcgctg
+
@4>C5?7FA?9@@@;1,-4B+E??4E9,DA>6CB)*2.)?7/.'65+75/
@h5n1_1
catcagttggattttgcctaaggatgtccaccatccttatcccaccaatt
+
7(9?)=C3AEF+@@F=C-E>EE57CD:;'?.*?7;>:3>2.;29..)=2'
@h5n1_2
catgcatttgttggaataatgttctcacaaatcaactgtattgaccgctg

This fixes the problem and allows BWA-MEM to work. However, header information is discarded. Is there another way to ensure the header names are the same while retaining all header information?

Assuming the read name lines only contain things like @JQ714243-920-1 and nothing else, then sed -i '1~4 s/-[12]$//' file.fastq (N.B., example updated a bit for safety, thanks to @terdon) will remove the -1 and -2 suffixes.

However, you should try to figure out why the read names were munged into that form to begin with (at least Illumina machines won't produce that, I don't know about others).


How can I dock a protein to a nucleic acid? I have a protein of interest and I would like to now how it interacts with RNA. I have structures of both molecules.

What tool can I use?

I'm not sure what you meant but you can take a look at NPDock (disclaimer we wrote that tool). If you have a structure of your protein of interest, you can dock it to the structure of your DNA/RNA of interest. Mind that this is a rigid body docking which means that the structure will not change upon binding.


Experimental Design for Differential expreression analysis I have a Normal esophageal Fibroblasts (NOFs) cultured in DMEM media; The same NOF also have been cultured with a tumor sample from a patient named 005 on DMEM media; I have also Cancer Associated Fibroblasts (CAFs) cultured on DMEM, on DMEM + tumor from a patient named as 005 and DMEM + tumor from a patient named 036 something like this picture;

I want to know if NOF and CAF are different in terms of gene expression but I am not sure about design for DESeq2 or Edger and ANOVA.

By @Devon's help I have done so

> mycols
   cell_type co_cultured
G2       NOF          NO
G3       NOF         YES
G4       CAF          NO
G5       CAF         YES
G6       CAF         YES

dds=DESeqDataSetFromMatrix(countData = NOFCAF,colData = mycols, design        =~ cell_type*co_cultured)

Now, I am wondering what would be a suitable contrast for results function in DESeq2

> resultsNames(dds)
[1] "Intercept"                   "cell_type_NOF_vs_CAF"        "co_cultured_YES_vs_NO"
[4] "cell_typeNOF.co_culturedYES"

> res=results(dds)
> res
log2 fold change (MLE): cell typeNOF.co culturedYES
Wald test p-value: cell typeNOF.co culturedYES
DataFrame with 2545 rows and 6 columns

Thanks a lot for correcting me

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

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


Integration of different microarray dataset to run GSEA I'm planning to run a GSE Analysis on some microarray datasets. Until now what I have done is to pre-processed them using frma algorithm.

The reason I chose frma instead of gcrma or rma is that I wanted to merge these different dataset into one.

Here is how I'm planing to proceed:

i.e

dataset1 : 3x Controls , 3x Treatment1
dataset2 : 2x Controls , 2x Treatment2
dataset3 : 2x Controls , 2x Treatment3

After merging

7x Controls , 3x Treatment1 , 2x Treatment2 , 2x Treatment3

Then I will execute the GSEA for controls , Treatmen1 , Treatment2 and Treatment3 separably.

So the question is if it's more reasonable to follow the path that I described or to run the GSEA on the separate datasets (for example for each control ) and then get the intersect of the gene sets of these three repetitions.

It makes more sense to evaluate by separate the pathway or gene set you want and see if in the three datasets result in a coherent message than to merge these datasets, as you will mix batch effects and treatment effects.

When in several datasets one looks for the same statistic it is usually considered meta analysis, not integration.


ATAC seq density calculation This paper Fig2b

b Changes of chromatin accessibility in AMD (n = 14) relative to
normal (n = 11) in all retina samples. Each dot represents one
ATAC-Seq peak. Blue line in the left panel indicates average fold
changes of peaks with the same ATAC-Seq intensity. The percentage of
reduced peaks is shown under the density curve in the right panel.

How is the density calculated or what is being shown there? average ATAC seq signal or Fold change?

They are not referring to density of reads. It is the density of the fold change, the y variable in the plot to left of Fig2b. It tells you the distribution of the foldchange and 92% of the fold change is <0.
See this page for more details on density and histograms
Suppose your data is something like this (purely making this up):
set.seed(555)
fc = rnorm(1000,-0.5,2)
averATAC = -0.5*(fc - 0.5)^2 + 0.2*fc + 3  + rnorm(100,25,10)

You fit a density to it and you get the plot:
par(mfrow=c(1,2))
plot(averATAC,fc)
plot(dens$y,dens$x,xlab = "density",ylab="fc")


Changing active.ident in Seurat Im trying to change the active.ident to another column in metadata but this error keeps popping up! I recently upgraded to R version 4.0.2 from 3.6.1 The older version was working but the new one isn't. I don't understand how to fix it?
M@active.ident <- M@meta.data$status
Error in (function (cl, name, valueClass)  :
assignment of an object of class “character” is not valid for @‘active.ident’ in an object of class
“Seurat”; is(value, "factor") is not TRUE

head(M@meta.data)
        orig.ident nCount_RNA nFeature_RNA percent.mito   status  RNA_snn_res.0.6 seurat_clusters
BC03_03       BC03   998414.7         3356   0.07869098 nonTumor                0               0
BC03_09       BC03   999853.6         7616   0.12964336  primary                0               0
BC03_13       BC03   998206.5         8457   0.17672205  primary                0               0
BC03_17       BC03   999535.1         7494   0.11323008  primary                0               0

I want to run the UMAP with the clusters named according to the status and not seurat clusters.However, the code isnt working..

M <- SetIdent(M, value = "status")
or more explicitly
M <- SetIdent(M, value = M@meta.data$status)
You can also use the group.by argument of UMAPPlot() or other plotting functions from Seurat for that matter.


How can I count the number of reads that support a variant in a bam file? I am calling variants from a human sample using bwa mem to align the reads and gatk to call the variants. I'm trying to understand why a specific variant was not called in my sample. I have checked the bam alignments in a GUI viewer and I can see that there are reads supporting the missing variant. It looks like the issue is a low allelic balance, with far more reads supporting the reference than the alternate allele but I want to get the actual numbers.

So, given a specific variant like this:

chr22   425236  C T

How can I count the number of reads in my sample.bam file that support that variant and the number that don't on Linux?

If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically:

samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam |
    cut -f 5 | tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c

That'll give you a count of how many of each base were seen.


Reconstruct a protein's 3D structure from a pair-wise distance matrix What is a good method to reconstruct a protein's 3D structure from a pair-wise distance matrix of residues?

I think the question is self-explanatory. I know of this method
but it seems that it uses contact maps, which actually use a certain cut-off. I don't have a cutoff based contact matrix but a matrix with actual distances.

Finding coordinates of N points that preserve (as well as possible) pair-wise distances from an NxN matrix is called multidimensional scaling.

But I don't enough about the context of your problem to tell if/how it is applicable.


How to run trimmomatic in HPC Sorry if this has been asked before, I've done a quick search and I don't think there's an easy explanation for me to understand.

I'm really new to bioinformatics/RNASEQ analysis, having only taken up an attempt to process my RNASEQ data last week. I have 0 background and some rudimentary knowledge of command line. So far I've managed to get access to my universities HPC and perform a fastqc on my fastq data.

What I'm stuck at at the moment is attempting to trim my RNASEQ data and remove adapters. My data had been run through Novaseq and I'm not exactly sure how to trim my current RNASEQ data.

I've attempted to run the following script using trimmomatic

# go to working folder
cd /project/CHR_RNASEQ/reads

# run trimmomatic
trimmomatic SE -threads 8 -phred33 SLIDINGWINDOW:4:20 ILLUMINACLIP:$TRIM/adapters/TruSeq3-SE-2.fa:2:30:10 MINLEN:36 1_R1.fastq.gz 1_R1_trimmed.fastq.gz

It did log a job in the HPC, but I had no output at the end. I didn't receive an error either. This is unfortunately just some code I had lifted and trialed and errored until something finally ran, but I'm very sure I'm very wrong and any help with understanding this would be really appreciated.

Cheers

Depending on the cluster management tool, you might have received e-mails when the "job" begins and ends. If so, you can check the "Exit status" of the job. For example, in the case of our HPC the relevant lines from the "job completion e-mail" would be:

Execution terminated
Exit_status=0

BTW, at first I was confused with the phrase "Execution terminated", however, in the scheduler used in our system, it means "the run has been successfully or unsuccessfully completed". Exit_status=0 int the next line means a successful run and for the unsuccessful runs there is a list with the error codes.

Secondly, if your run was successful, you would see new fastq files containing "_trimmed" in the filename (reading this from your code). These fastqs should be more or less similar in size to the original depending on how aggressive the trimming settings are.

If you can provide more information, i.e. the type of cluster management tool used, the log file, ..., you can get more help and precise answers from the community.


Basic question about finding mRNA sequence in transcriptome I am a computer scientist just starting out in bioinformatics topics, and would appreciate any guidance that can be given here:
I have an mRNA sequence- an isoform - whose length is about 4000 base pairs (bp) - the origin is provided in the link.
I also have a bam file/fasta file of the entire transcriptome of a kidney cell (see Data Access tab).
My objective is to quantify mRNA expression levels of the mRNA sequence in the transcriptome.
How does one get started? Should I use the original format (TenX bam, about 10 GB) or SRA archive data (fasta, about 3.5 GB)?  Why are the files so different in size?
I have a MacBook with samtools installed and 16GB of RAM. Is this enough to process the data and what basic commands should be used to quantify the 4000 bp sequence overlaps in the transcriptome?

Basically you want to know the VHL Gene expression in that 10x dataset. They have provided the barcode, matrix and feature files on GEO GSE131685 as Supplementary file. You could use those files and get the gene expression matrix as described on 10x website.
library(Matrix)

matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")

mat <- readMM(file = matrix.path)

feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)

barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

This mat can be further loaded into Seurat for any downstream analysis. All this can be done on a laptop.
If you want to re-quantify the gene expression for whatever the reason, you need to re-run the cellRanger pipeline using GTF of your interest (it needs some preparation), which may not work on your laptop. It requires large memory and long running times, ideally to submit as job on cluster.


levels of factors in the design have non-unique level names after make.names() is applied When I try to import the data using DESeq2 package:

dds <- DESeqDataSetFromTximport(txi, sampleTable, ~Group + Cre)

I am getting the error:

  Error in validObject(.Object) :
    invalid class “DESeqDataSet” object: levels of factors in the design have non-unique level names after make.names() is applied.
    best to only uobject letters and numbers for levels of factors in the design

Here is how my sampleTable looks like:

    ID      Group Cre
1  Bla_36   IL4   -
2  Bla_38   IL4   -
3  Bla_11   IL4   -
4  Bla_17   IL4   -
5  Bla_20   IL4   +
6  Bla_21   IL4   +
7   Bla_8   IL4   +
8  Bla_22   IL4   +
9   Bla_7   IL4   +
10 Bla_35   PBS   -
11 Bla_37   PBS   -
12 Bla_13   PBS   -
13  Bla_1   PBS   -
14  Bla_9   PBS   +
15 Bla_23   PBS   +
16 Bla_19   PBS   +
17 Bla_24   PBS   +

What am I doing wrong here? And how to make it working?

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

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

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

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


Extracting concatenated metadata from a column in R I've downloaded some ChIP-seq data from www.chip-atlas.org and I want to use this to conduct some enrichment tests with genomic regions from my experiment. The .bed file generated by that website is set up so that each row represents a genomic region from a specific ChIP-seq experiment. I have the data set up as a dataframe in R.
Each row has a single column with all of the metadata, including such important information as the target transcription factor, model system, and tissue type.  These are delimited by semi-colons, but the number of pieces of information varies between rows so it's impossible to simply split by the delimiter.
The metadata is structured something like this:
name=CTCF; strain=B6; age=adult; tissue=brain
name=Mecp2; age=embryo; tissue=5t3 cells; passages=5-10

So within each row, the order and number of metadata terms varies.
I want to reshape this column into something like this:
name; strain; age; tissue; passages
CTCF; B6; adult; brain; NA
Mecp2; NA; embryo; 5t3 cells; 5-10

How can I extract and reshape this data in R? It seems like it should be possible, but the normal approaches to splitting columns by delimiters are failing me since the rows all have different numbers of items and order.
Here's a slice of the real data which is much more messy than my little example.
test_df <- data.frame(pos=c(3002983, 3002881, 3010946, 3021775), metadata=c("ID=SRX661585;Name=Ctcf%20(@%20Brain);Title=GSM1446329:%20CTCF%20brain%20ChIPSeq%3B%20Mus%20musculus%20x%20Mus%20spretus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Mouse%20adult%20hybrid%20brain;tissue=Adult%20whole%20brain;strain=BL6/spretus;chip%20antibody=CTCF,%20MillIpore%2007-729;", "ID=SRX2957352;Name=Smc1a%20(@%20AtT-20);Title=GSM2684855:%20SMC1%20ChIPseq%20in%20AtT-20%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=AtT-20%20cells;cell%20line=AtT-20%20cells;cell%20type=pituitary%20corticotroph%20cell;passages=10-20;chip%20antibody=Bethyl%20A300-055A;", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);"))

Thank you for your help!

In general you want to:

Clean your text
Split the features by semicolons
Split keys and values by '='
Create a temporary data frame with keys and value
Reshape

library(tidyverse)
library(urltools)
library(textclean)
test_df <- data.frame(pos=c(3002983, 3002881, 3010946, 3021775), metadata=c("ID=SRX661585;Name=Ctcf%20(@%20Brain);Title=GSM1446329:%20CTCF%20brain%20ChIPSeq%3B%20Mus%20musculus%20x%20Mus%20spretus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Mouse%20adult%20hybrid%20brain;tissue=Adult%20whole%20brain;strain=BL6/spretus;chip%20antibody=CTCF,%20MillIpore%2007-729;", "ID=SRX2957352;Name=Smc1a%20(@%20AtT-20);Title=GSM2684855:%20SMC1%20ChIPseq%20in%20AtT-20%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=AtT-20%20cells;cell%20line=AtT-20%20cells;cell%20type=pituitary%20corticotroph%20cell;passages=10-20;chip%20antibody=Bethyl%20A300-055A;", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);"))

tbl <- apply(test_df, 1, function(row){
  # Split by smicolons and clean html tags
  n <- replace_html(unlist(strsplit(row[2], ";")))
  # Decode urlencoded characters, mostly space
  n <- url_decode(n)
  # Split into keys and values
  s <- do.call(rbind, str_split(n, "="))
  # Get keys and values into vectors and clean up leading and trailing whitespace
  v <- s[, 2]
  v <- str_replace_all(v, "(^ +)|( +$)", "")
  k <- s[, 1]
  k <- str_replace_all(k, "(^ +)|( +$)", "")
  # Create final data frame with original columns (pos), and also columns from
  # keys and values
  tibble(pos = row[1], key = k, value = v)
})

# Concatenate list into tibble
(tbl <- do.call(rbind, tbl))

This tibble looks like:
# A tibble: 37 x 3
   pos     key           value
   <chr>   <chr>         <chr>
 1 3002983 ID            SRX661585
 2 3002983 Name          Ctcf (@ Brain)
 3 3002983 Title         GSM1446329: CTCF brain ChIPSeq; Mus musculus x Mus spretus; ChIP-Seq
 4 3002983 Cell group    Neural
 5 3002983 source_name   Mouse adult hybrid brain
 6 3002983 tissue        Adult whole brain
 7 3002983 strain        BL6/spretus
 8 3002983 chip antibody CTCF, MillIpore 07-729
 9 3002881 ID            SRX2957352
10 3002881 Name          Smc1a (@ AtT-20)
# … with 27 more rows

Using this, you can:
# Pivot wider
tbl %>%
  pivot_wider(names_from = key, values_from = value)

Which nets you:
# A tibble: 4 x 15
  pos   ID    Name  Title `Cell group` source_name tissue strain `chip antibody` `cell line` `cell type`
  <chr> <chr> <chr> <chr> <chr>        <chr>       <chr>  <chr>  <chr>           <chr>       <chr>
1 3002… SRX6… Ctcf… GSM1… Neural       Mouse adul… Adult… BL6/s… CTCF, MillIpor… NA          NA
2 3002… SRX2… Smc1… GSM2… Neural       AtT-20 cel… NA     NA     Bethyl A300-05… AtT-20 cel… pituitary …
3 3010… SRX5… Mecp… GSM3… Neural       Olfactory … NA     NA     MeCP2 antibody… NA          Olfactory …
4 3021… SRX5… Mecp… GSM3… Neural       Olfactory … NA     NA     MeCP2 antibody… NA          Olfactory …
# … with 4 more variables: passages <chr>, `strain/background` <chr>, `genotype/variation` <chr>, age <chr>


How to find adapter sequence I have this GSE dataset ( GSE104279 ) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104279)
I want to run cutadapt, but how would I find the adapter sequence so I can run cutadapt?

The definitive answer is the sequencing institute/lab. They know what protocol/chemistry etc. was used.
If you don't have access to that a number of tools check for known adapter sequences. Run e.g. FASTQC, which will tell you the proportions of adapter sequence found. Tools like TrimGalore can also autodetect the most common adapters.
UPDATE:
I checked 2M reads from one run of your dataset and it seems to me that it was uploaded already preprocessed:

Check all your samples, but I don't see a need to do any preprocessing if all of them look like this.


Is there a safe catch-all adapter sequence for trimming? I would like to trim/mark adapters using trimmomatic or picard MarkIlluminaAdapters from a series of Illumina Paired-End read fastqs. The fastq files may have been done using different kits or different generations/versions of illumina chips.

Without knowing ahead of time how the libraries were prepared (or which kit was used), Is there a safe "catch all" adapter sequence (or set of such sequences) that can be passed to a trimming tool and that would cover all Illumina kits? In other words, what are we expected to know about the fastq to ensure correct trimming?

edit: Picard defines a few built-in adapters which cover many standard kits, but it's not immediately obvious if you need to pass a common one (built-in sequences defined here).

You're best off just using fastp or Trim Galore!, both of which will determine the adapter sequence for you. Trim Galore! uses a built-in list of known sequences for this, whereas fastp uses read overlap.


How to read structural variant VCF? The IGSR has a sample for encoding structural variants in the VCF 4.0 format.

An example from the site (the first record):

#CHROM  POS   ID  REF ALT   QUAL  FILTER  INFO  FORMAT  NA00001
1 2827693   . CCGTGGATGCGGGGACCCGCATCCCCTCTCCCTTCACAGCTGAGTGACCCACATCCCCTCTCCCCTCGCA  C . PASS  SVTYPE=DEL;END=2827680;BKPTID=Pindel_LCS_D1099159;HOMLEN=1;HOMSEQ=C;SVLEN=-66 GT:GQ 1/1:13.9

How to read it? From what I can see:

This is a deletion (SVTYPE=DEL)
The end position of the variant comes before the starting position (reverse strand?)
The reference starts from 2827693 to 2827680 (13 bases on the reverse strand)
The difference between reference and alternative is 66 bases (SVLEN=-66)

This doesn't sound right to me. For instance, I don't see where exactly the deletion starts. The SVLEN field says 66 bases deleted, but where? 2827693 to 2827680 only has 13 bases between.

Q: How to read the deletion correctly from this structural VCF record? Where is the missing 66-13=53 bases?

I just received a reply from 1000Genomes regarding this. I'll post it in its entirety below:

  Looking at the example you mention, I find it difficult to come up with an
  interpretation of the information whereby the stated end seems to be correct,
  so believe that this may indeed be an error.

  Since the v4.0 was created, however, new versions of VCF have been introduced,
  improving and correcting the specification. The current version is v4.3
  (http://samtools.github.io/hts-specs/). I believe the first record shown on
  page 11 provides an accurate example of this type of deletion.

  I will update the web page to include this information.

So we can take this as official confirmation that we were all correct in suspecting the example was just wrong.


HMM Profile from Convergence of retrotransposons in oomycetes and plants I am trying to recreate the experiment from Convergence of retrotransposons in
oomycetes and plants by Kirill Ustyantsev, Alexandr Blinov and Georgy Smyshlyaev. In the "Methods" section a

  HMM profile was constructed from the amino acid alignment of aRNH
  sequences from the Ustyantsev et al

I was wondering where I could get this profile. I checked the cited papers but they were also referencing other papers. Could someone maybe clarify where I could get the profile?

The alignment they used to build the HMM is included as supplementary file 8 in the paper you cite:

  Multiple amino acid sequence alignment of RNH genes and domains from diverse taxa constructed and used for the phylogenetic reconstruction in the present study. (TXT 45 kb)

So, you can download that alignment and then use a tool like hmmbuild (part of HMMER) to build the HMM profile:

hmmbuild model.hmm 13100_2017_87_MOESM8_ESM.txt

The command above will create the file model.hmm which contains the HMM built from the alignment in the supplementary file 13100_2017_87_MOESM8_ESM.txt. I am including the first few lines here (this is not the entire file, don't use this as the model—build it yourself. I only include it to illustrate the format):

HMMER3/f [3.1b2 | February 2015]
NAME  13100_2017_87_MOESM8_ESM
LENG  121
ALPH  amino
RF    no
MM    no
CONS  yes
CS    no
MAP   yes
DATE  Mon Dec 11 13:33:29 2017
NSEQ  165
EFFN  5.697556
CKSUM 2893352190
STATS LOCAL MSV       -9.8975  0.71332
STATS LOCAL VITERBI  -10.6450  0.71332
STATS LOCAL FORWARD   -3.9247  0.71332
HMM          A        C        D        E        F        G        H        I        K        L        M        N        P        Q        R        S        T        V        W        Y
            m->m     m->i     m->d     i->m     i->i     d->m     d->d
  COMPO   2.34473  4.21240  2.85435  2.69203  3.54228  2.79153  3.54812  3.01140  2.72630  2.48370  3.82225  2.94269  3.49156  3.05882  2.84632  2.63343  2.86306  2.72633  4.27012  3.56181
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.13773  5.45955  2.08419  0.61958  0.77255  0.00000        *
      1   2.99488  4.57852  4.44659  3.29307  1.80597  4.16203  3.88788  2.87265  3.75733  2.38594  3.30344  3.33633  4.53284  3.64324  3.69630  2.90350  3.30742  2.83413  4.36991  1.50673      1 y - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00723  5.32905  6.05140  0.61958  0.77255  0.45275  1.01027
      2   2.66602  2.17750  5.01715  4.40788  1.72617  4.29647  4.62737  2.89579  4.21116  2.71668  3.64495  4.43195  4.66188  4.35829  3.99238  3.10363  2.02207  1.89528  5.14692  3.02155      2 f - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00661  5.41881  6.14115  0.61958  0.77255  0.61842  0.77391


Control width of samtools tview "snapshot" when redirecting I have a large number of loci I would like to examine manually with samtools tview. Rather than typing or copy-n-pasting dozens of coordinates, I was hoping to automate the generation of pileup "images" (text files).

With some experimenting just now, I was happy to learn that if I specify the display option -d H or -d T the output can easily be redirected to a file for later inspection. However, it appears that I cannot control the window captured. Changing the width of the terminal window doesn't change the output, nor does adjusting the -p flag (i.e. -p chr1:501-700). No matter what, tview always reports an 80bp window in text or HTML mode.

Is there any way I can adjust this?

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

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


Converting hmmer --tblout to pandas dataframe Is there a way to convert a hmmer output to a pandas dataframe? I am also unsure how to load a hmmer tblout table into python via the Bio module. I believe you can call a hmmer format with SeqIO.parse or SeqIO.search. How can you convert the table to a pandas dataframe?

The manual says fields in a tabular output file are space-delimited. It should be fairly straightforward to load a space-delimited file in pandas.

From the manual:

  --tblout <f> Save a simple tabular (space-delimited) file summarizing the per-target output, with one data line per homologous target model found.


How to reverse complement the DNA sequences for given inverse/reverse coordinates? I have the series of coordinates in id.txt file, whose coordinates sequences are in genome.fasta file.
The coordinates of id.txt file is shown below,

Contig3:15-7
Contig2:5-10
Contig1:12-3

The genome.fasta file is shown below,

>Contig1
AAGGCCATCAAGGACGTGGATGAGGTCGTCAAG
>Contig2
ACGAAGCGCTCGCCAAGGCCGAAGAAGAAGGCC
>Contig3
GCTGCGGCGCTGATCCTGGCGGCCCGCGCCGAG

I have used the following codes to extract the sequences from genome.fasta based on the coordinates in id.txt by using the following command xargs samtools faidx genome.fasta <id.txt > result.fasta
It has extracted only Contig2:5-10sequence. Because, the coordinates of Contig2:5-10is in proper order. But Contig3:15-7 and Contig1:12-3 coordinates are in inverse/reverse order, So samtools could not fetch those sequences.

I need to extract those sequences from Contig1:12-3 and Contig3:15-7 coordinates and also I want to reverse complement them. I have  larger fasta file and plenty of coordinates sequences to be extracted, So I need to automate this process. So, please help me to automate this process.

The expected output is

>Contig3
GATCAGCGC
>Contig2
AGCGCT
>Contig1
CTTGATGGCC

Thank you in advance.

Since you tagged python, here is a solution that solves this in a single script using Biopython:

import sys
from Bio import SeqIO

d = {}
with open('id.txt') as f:
    for line in f:
        contig, coords = line.split(':')
        start, end = map(int, coords.split('-'))
        d[contig] = start, end

for record in SeqIO.parse('genome.fasta', 'fasta'):
    start, end = d[record.id]
    if start > end:
        record.seq = record.seq[end - 1: start].reverse_complement()
    else:
        record.seq = record.seq[start - 1: end]
    SeqIO.write(record, sys.stdout, 'fasta')

Output on your example:

>Contig1
CTTGATGGCC
>Contig2
AGCGCT
>Contig3
GATCAGCGC


How do you write a .gz fastq file with Biopython? How do you write a .gz (or .bgz) fastq file using Biopython?

I'd rather avoid a separate system call.

The typical way to write an ASCII .fastq is done as follows:

for record in SeqIO.parse(fasta, "fasta"):
    SeqIO.write(record, fastq, "fastq")

The record is a SeqRecord object, fastq is the file handle, and "fastq" is the requested file format. The file format may be fastq, fasta, etc., but I do not see an option for .gz.

Here is the SeqIO API.

I'm not sure I'm doing it the best way, but here is an example where I read a compressed gzip fastq file and write the records in block gzip fastq:

from Bio import SeqIO, bgzf
# Used to convert the fastq stream into a file handle
from io import StringIO
from gzip import open as gzopen

records = SeqIO.parse(
    # There is actually simpler (thanks @peterjc)
    # StringIO(gzopen("random_10.fastq.gz").read().decode("utf-8")),
    gzopen("random_10.fastq.gz", "rt"),
    format="fastq")

with bgzf.BgzfWriter("test.fastq.bgz", "wb") as outgz:
    SeqIO.write(sequences=records, handle=outgz, format="fastq")