Fix dash/empty ALT alleles for deletions in VCF I have a VCF file that has deletions specified like this, with the REF allele containing only the deleted portion and the ALT allele just saying -:

chr1 101 . AAA -

I would like to combine this with another VCF in which deletions are specified like this, with REF and ALT both extended by 1bp to the left to ensure neither is empty:

chr1 100 . GAAA G

Is there a tool to convert from one of these representations to the other? I have already tried vt normalize and bcftools norm but neither of them seem to do this.

I don't really mind which representation I end up with although I think the second is the one more consistent with the spec.

According to the specs the ALT is not allowed to contain -.
Here is a little python (>=3.6) script that is able to fix it. It makes use of the pysam module. You need a reference file and this fasta file must be indexed by samtools faidx reference.fa. Please note that it doesn't support multiple values in the ALT column for records that contain a -.
Usage:
python fixvcf.py -f reference.fasta -o out.vcf in.vcf

The script:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import sys

import pysam

class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """
        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) -> pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """
        if "-" in record.alts:
            ref = self.fasta.fetch(region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1})")
            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "-" else f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see: https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """
        print(self.vcf.header, end="")
        print(*map(self.fix_record, self.vcf), sep="", end="")

if __name__ == "__main__":
    app = Application()
    app.start()

A lot of the code above is boilerplate for running a script on the command line. The critical code for the transformation is in the fix_record function. The relevant code is below, expanded, with line numbers and a description of each command.
[1]   ref = self.fasta.fetch(
          region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1})"
      )
[2]   record.pos = record.pos - 1
[3]   record.ref = f"{ref}{record.ref}"
[4]   record.alts = tuple(
          map(
              lambda alt: ref if alt == "-" else f"{ref}{alt}",
              record.alts
          )
      )

This grabs the base preceding the inserted/deleted sequence.
This corrects the position of the variant to account for prepending another nucleotide.
This defines the new corrected REF allele: concatenate the base grabbed in [1] with the original REF allele
Apply an anonymous (lambda) function to each ALT allele; the function will replace - with the base grabbed in [1], or preprend the base to the allele if it is not -


Test to determine if two genes/exons share the same evolutionary histories? In classic phylogenetic inference one is usually given various orthologue sequences of a given gene across various species. Those sequences are then multiple aligned and used to construct a phylogenetic tree.

Say, I split the multiple alignment of this gene into two parts (part A and B).
My question now is: What is the most appropriate test to determine if part A and B share the same phylogenetic history?

Having looked through literature, it appears this can be done by constructing phylogenetic trees separately for part A and B and test for incongruence using

Kishino–Hasegawa (KH test);
Shimodaira–Hasegawa (SH-test);
SOWH-test etc.

as nicely summarized in Plantet et all, and Goldman et al. etc.

However the above test rely on constructed trees beforehand which might be incorrect due to inappropriate algorithms/parameters.

I was therefore wondering if there are more appropriate test maybe without having to construct trees at all?

Since there hasn't been any other answers, I'd thought to share my approach in case someone is having the same issue.

Gauging from the comments of tendon and NatWH, it seems using the SH-test valid stategy to tackle this problem.

Having done some research, I could only find one functional software package, the R-package phangorn which implements both, the SH- as well as SOWH-test and thought to share my R-workflow:

Reading in multiple-alignment (MLN). Here: amino-acids (AA)

 library(phangorn)
 mln <- read.phyDat("my_mln.fasta", format = "fasta", type = "AA")

Partitioning MLN into two parts (part A and B) given a predefined breakpoint (here 100)

breakpoint <- 100
 # partA: 1:breakpoing
 ( mln_partA <- subset(mln, select = 1:breakpoint, site.pattern = FALSE) )
 # 58 sequences with 100 character and 45 different site patterns.
 # The states are a r n d c q e g h i l k m f p s t w y v

 ## partB: breakpoint:end
 ( mln_partB <- subset(mln, select = (breakpoint + 1):length(attr(mln, "index")), site.pattern = FALSE) )
 # 58 sequences with 1702 character and 720 different site patterns.

Constructing tree(s) separately from part A and B using Neighbourhood Joining (NJ) method

 ## Note, other tree building algorithms might be used
 ( treeNJ_partA <- NJ(dist.ml(mln_partA)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.
 ( treeNJ_partB <- NJ(dist.ml(mln_partB)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.

Using topology of previously constructed trees and compute optimal likelihood based on partA data

 ( treeNJ_partA_fit1 <- optim.pml(pml(treeNJ_partA, mln_partA)) )
 # optimize edge weights:  -1402.126 --> -1379.043
 #
 #  loglikelihood: -1379.041
 #
 # unconstrained loglikelihood: -302.6292
 ( treeNJ_partB_fit1 <- optim.pml(pml(treeNJ_partB, mln_partA)) )
 # optimize edge weights:  -1633.646 --> -1530.745
 #
 #  loglikelihood: -1530.682
 #
 # unconstrained loglikelihood: -302.6292

Performing SH-test for for both trees

 ( shPart1_t12 <- SH.test(treeNJ_partA_fit1, treeNJ_partB_fit1, data = mln_partA, B = 1e5) )
 #      Trees      ln L Diff ln L p-value
 # [1,]     1 -1379.041     0.000 0.48461
 # [2,]     2 -1530.682   151.641 0.00000

As can be seen in this example, the tree constructed for partB doesn't fit partA data as much as the tree constructed from partA.


Change sample names while using cummeRbund I want to analyze RNA seq data using R package cummeRbund. Using this package I have a problem changing sample names. When I run this command in cummeRbund

samples(cuff)
sample_index sample_name sample_name parameter value
1            1          q1        <NA>      <NA>  <NA>
2            2          q2        <NA>      <NA>  <NA>
3            3          q3        <NA>      <NA>  <NA>
4            4          q4        <NA>      <NA>  <NA>
5            5          q5        <NA>      <NA>  <NA>

It shows q1,q2,q3,q4 and q5 as sample names for the five samples but I want to change the sample names and put following names for the samples Batis, Blue_silver, Chakwal_50, Local_white and Syn_22 for q1,q2,q3,q4 ad q5 respectively. To accomplish this I ran following command:

> samples(cuff) <- c("Batis", "Blue_silver","Chakwal_50","Local_white","Syn_22")
Error in samples(cuff) <- c("Batis", "Blue_silver", "Chakwal_50", "Local_white",  :
could not find function "samples<-"

but this command gave error as shown above.

Also the BAM files have same names for all the samples as mentioned above.

How could I change the sample names to desired names in R using cummeRbund?

Or I need to re-run cuffdiff and change the names somewhere else?

Don't know which file was used by cummeRbund for sample names but replacing q1,q2,q3 and q4  with Batis, Blue_silver, Chakwal_50 and Syn_22 for all the files in the directory containing cuffdiff output worked for me.

I used sed language for this as shown below:

sed -i 's/q4/Syn_22/g' *

Or more generally you can use

sed -i 's/sample_name_shown/sample_name_to_change/g' *

After that I deleted the sqlite object created by cummeRbund and again created this object using readCufflinks() in R.
This time sample names are changed as desired:

samples(cuff)
sample_index sample_name sample_name parameter value
1            1       Batis        <NA>      <NA>  <NA>
2            2 Blue_silver        <NA>      <NA>  <NA>
3            3  Chakwal_50        <NA>      <NA>  <NA>
4            4      Syn_22        <NA>      <NA>  <NA>


awk 'BEGIN{OFS="\t"}{if($1==lchrom && $4==lname && $2 == lend) {lend = $3}else{if(lchrom) {print lchrom, lstart, lend, lname;}; lchrom=$1; lstart=$2; lend=$3; lname=$4}}END{print lchrom, lstart, lend, lname}' foo.bed > new.bed

That produces:

chr1    0   10000   LOC101929784::chr1:13201-13800::E1
chr1    10000   10600   LOC101929784::chr1:13201-13800::XYZ
chr1    10600   10800   LOC101929784::chr1:13201-13800::E1
chr1    10800   13200   DDX11L1::chr1:11868-14409.E1
chr1    13200   29600   WASH7P::chr1:14403-29570.E1
chr1    29600   31400   MIR1302-2HG::chr1:29553-31109.E1

This explicitly checks for bookends and nothing else and assumes that everything is sorted.


What's "signal" in gene expression? I'm exploring gene expression on public available databases. RefExa is one of the databases I've been studying.

When I searched something there, the website reported the expression value unit as "signal". See screenshot where I highlighted:

I'm not able to find anywhere what that really mean. TPM? RPKM? Raw counts? Normalizated for gene length?

My question is that how to interpret gene expression "signal"? What does signal 100 mean?

Looking at this resource, these data are microarray data rather than RNAseq data. Thus the "signal" in question is the intensity of florescent light being emitted from the the probe associated with your gene.  The relationship between intensity and quantity of RNA is different for every probe because each probe will have a different affinity for complementary fragments dependent on the sequence in question (not least the melting temperature).

Affy intensities are usually presented as log values going up to around 12, but the raw number is next to meaningless and usually only means something as a normalized (e.g. RMA) log fold change.

I've never come across RefExA before, and it appears to be undocumented and unmaintained. I wouldn't rely on it for usable data if I were you. A more modern database of expression data is the EBI Expression Atlas. It presents data from a wide range of different studies, processed and normalised on a unified pipeline and expressed as TPMs.


What is a good pipeline for using public domain exomes as controls? I'm currently attempting association analysis with an extremely small set of patient exomes (n=10), with no control or parental exomes available. Downloading the ExAC VCF of variant sites (http://exac.broadinstitute.org/downloads) or the 1000G integrated call sets (http://ftp.1000genomes.ebi.ac.uk/) and combining this with our pooled patient VCFs has not been successful (I suspect the approach of attempting to merge such large VCFs generated from different pipelines is rather naive).

Looking at the primary literature, I have gathered it should be possible to use these resources to help increase statistical power for our analysis. My question is how do I take these large .vcfs with many samples and successfully merge them to our patient .vcfs, such that the combined VCF can be used downstream to run analysis packages? (PODKAT, PLINK, etc.)

I gather that you want to use the background SNP frequency as a prior for input to your SNP calling algorithm?

I'm not sure of a canned algorithm for doing this, but a quick google shows up some promising links:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3593722/
https://software.broadinstitute.org/gatk/documentation/article.php?id=4723
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0489-0


Extracting modules from correlation and then module clustering This is the paper where they have first done co-expression analysis of TF and then found coexpressed modules which were further subjected to module clustering ,paper

So now im clear about the first part of doing the coexpression ,the part which they have written on the paper Im quoting

  Using Pearson correlation coefficient of 0.4 as a cutoff to define a
  positive correlation between TF pairs, we obtained 37 tightly
  related TF modules that showed excellent co-expression patterns (Fig.
  4a and Supplementary Data 4).We summed the average intensity of each
  module and then performed hierarchical clustering to reveal the
  relationship among the 37 TF modules

As of now If i have done something like this ,i use all the modules coming up with no threshold as such which they have written 0.4, if i get it, then the first figure is a gene-gene pair correlation, how to filter it ,because if i get try to set a threshold i will have to look for every pair irrespective of the module found ,

How to do it in R im not sure about it, any suggestion would be really helpful,next one is how are they going for the next one which is module clustering ,as there are 37 TF modules ,each of them having different gene distribution ,here the second figure the modules are in the rows if i understand ,so how have they done it i would like to get an idea about it.

Earlier i have asked a question , i have the data as if that can be used to show me how it can be done that would be really helpful.

Co-expression in R is doable via gglot2 at the very least... I see if I can find the code.

Essentially you are performing a co-variance analysis, however Pearson' at 0.4 is a LOW cut-off and it should be 0.7.

I've got a lot on right now, so it will take a week or so before I foward the relevant code. R will have other methods outside ggplot2 however.

I was thinking ggcorr() and ggpairs() which both in the GGally package of R-cran. The output is very cute, well its ggplot2 which is trendy.

library (ggplot2)
library (GGally)
library (dplyr)

df <- select(mydata, 1:5)
cor(df) # raw output
ggcorr(df)

output: heatmap similar to Fig. a

ggpairs(df)

output: very fancy matrix plot Its pairwise scatterplots, diagrams and numbers, which is not represented in Fig. a to c, but its a very fancy version of Fig. d ... Its simplistic analysis, but its a very cute plot.

ggcorr uses pairwise (obviously) Pearson by default, which is the most versatile.

I can definitely say that 0.7 for Pearson for a straight correlation (not a covariance analysis) is conservative low and as a reviewer I would immediately complain about anything lower. A heatmap is different depending on its interpretation.

The multivariate analysis of Pearson pairwise correlation, which is what you are requesting is a different question. UPGMA style analysis directly on the heatmap, is fine if there is "one heatspot", but otherwise has limitations IMO.


split 1 column input into 5 column bed file In f1 below I am trying to split $1 based on each line and create a bed file. If the line is a snp then the : is spilt and the text is $1 the last digit is $2 -1 and $3 and the letter to the left of the > is $4 and the letter to the right is $5. (lines 1 and 2)

If the line has dup then the : is spilt and the text is $1 the last digit is $2 and $3 and the $4 is - and the letter to the right is $5. (line 4)

If the line has del with nothing else then the : is spilt and the text is $1 the last digit is $2 and $3 and the letter to the left of the > is $4 and $5 is -. (line 5)

I am not sure how to format for line 3 (the complex indel)... maybe split $1 on the : and the text is $1 the last digit is $2 and $3 and the letters to the left of the > is $4 and the letters to the right is $5. Thank you :)

f1

chr7:140453145A>T
chr7:140453136A>T
chr7:140453135_140453136delCAinsTT
chr20:31022287dupA
chr19:13054614delG

desired

chr7    140453144   140453145   A   T
chr7    140453135   140453136   A   T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -

awk

awk 'BEGIN{OFS="\t"} {sub(/[^0-9]+$/, "", $1)}' f1

awk 'BEGIN { FS = "[ -]"; OFS="\t" } NF==3 { print $1, $2 - 1, $2, $3 } NF==4 { print $1, $2, $3, $4 }' f1

The following Python script seems to do the job.

#!/usr/bin/env python3

from __future__ import print_function
import re
import sys

for line in sys.stdin:
    if '_' in line:
        match = re.search(r'(\S+):(\d+)_(\d+)del([ACGT]+)ins([ACGT]+)', line)
        print(match.group(1), match.group(2), match.group(3), match.group(5), match.group(4), sep='\t')
    elif 'dup' in line or 'del' in line:
        match = re.search(r'(\S+):(\d+)(del|dup)([ACGT]+)', line)
        if match.group(3) == 'del':
            print(match.group(1), match.group(2), match.group(2), match.group(4), '-', sep='\t')
        else:
            print(match.group(1), match.group(2), match.group(2), '-', match.group(4), sep='\t')
    else:
        match = re.search(r'(\S+):(\d+)([ACGT]+)>([ACGT]+)', line)
        coord = int(match.group(2))
        print(match.group(1), coord - 1, coord, match.group(3), match.group(4), sep='\t')

Invoked like so on the command line.

[standage@lappy ~] $ ./transform < f1
chr7 140453144 140453145 A T
chr7 140453135 140453136 A T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -


DESeq2 multiple treatments, multiple time points, multiple cell lines Note: this question has also been asked on Bioconductor Support
I know there are a lot of questions asking similar things here on this forum and I checked the vingette and did a lot of other research but I still have a hard time wrapping my head around it.
In my position, data is unfortunately just put in front of me and I'm asked for pretty pictures. Even making it clear to the wetlab people in the lab that we need replicates is an issue.
I have multiple cell lines, multiple time points and multiple treatments:
Cell lines: CL1, CL2, CL3 Time points: 6h, 24h Treatments: T1, T2, T3, Control
For each cell line and each time point, there are 3 different treatments plus a control. 3 replicates for each sample -> 72 samples
What I want is actually quite simple; I want to measure control vs each treatment at each time point in each cell line. No testing across time points or across cell lines.
My first thought was to separate the data into 6 different data sets (CL1, 6h | CL1, 24h | CL2, 6h | CL2 24h | CL3, 6h | CL3, 24h) and simply do the DE analysis separately, but I read that this is not the way to go. In the end it would also be nice to get a normalized count matrix with all data normalized together, for PCAs and similar.
I hope I made it understandable. How do I design the DESeqDataSet?

To compare a subset of samples to another subset of samples, make one column of column data which will distinguish them all the way you want (so probably join cell line, treatment, timepoint), and make that column your design, and then when you call results, you specify the pair you want to contrast, one pair at a time.


Efficiently concatenating disjoint, ordered bigBed files Given a set of bigBed files, each of which are disjoint to each other and in lexicographical ("sort-bed") sort order, I would like to know if there is a way to concatenate these archives efficiently.

Presuming I know the order in which to specify files (to constrain the order of intervals), is there a way to concatenate them directly, i.e., without first extracting them to bed, concatenating the extracted files, and remaking a bigBed file?

No, there's no simple way to do that without someone writing a program tailored toward it that handles the footer section and recreates the zoom levels that follow the compressed blocks of intervals. If you wanted to write such a program, modifying bigWigCat from Kent's tools would be a reasonable tack.


Filtering step for read counts data I have around 1200 samples as columns and 60,000 genes with Htseq-Counts data. Before normalization with voom function I want to do filtering step.

I want to remove genes whose expression is == 0 in at least 10 samples.

Can I do this with read counts itself or should I transform counts to cpm?

If filtering with counts data is fine May I know how to do this.

If you have your counts in a data.frame called counts, something like this might work:

filtered.counts <- counts[rowSums(counts==0)<10, ]

For example lets assume the following data frame.

> A <- c(0,0,0,0,0)
> B <- c(0,1,0,0,0)
> C <- c(0,2,0,2,0)
> D <- c(0,5,1,1,2)
>
> counts <- data.frame(A=A, B=B, C=C, D=D)
> counts
  A B C D
1 0 0 0 0
2 0 1 2 5
3 0 0 0 1
4 0 0 2 1
5 0 0 0 2

And now you would want to remove the genes with at least 3 zeroes or more:

> filtered.counts <- counts[rowSums(counts==0)<3, ]
>
> filtered.counts
  A B C D
2 0 1 2 5
4 0 0 2 1

Every gene with 3 or more zeroes is now removed.


HIVdb program and FASTA sequence input I am trying to use HIVdb(https://hivdb.stanford.edu/hivdb/by-sequences/) program. I have a HIV protease which I want to test with most commonly used PIs. But when I entered it in fasta sequence format to the HIVdb, it give me error code saying

  Critical: There were no Protease, Reverse Transcriptase, or Integrase genes found, refuse to process.

For example, if I download FASTA sequence of 1HPV from PDB and paste it in HIVdb input sequence window as below,

PQITLWQRPLVTIKIGGQLKEALLDTGADDTVLEEMSLPGRWKPKMIGGIGGFIKVRQYDQILIEICGHKAIGTVLVGPT
PVNIIGRNLLTQIGCTLNF

I get the error code. Any suggestion on how to input sequences in HIVdb?
(Note: I am new to this subject)

The HIVdb page states the following in the second paragraph (emphasis mine):

  Sequences can be entered as plain text if just one sequence is entered. Sequences must be entered using the FASTA format if multiple sequences are entered. Sequences can be pasted in the text box or uploaded using the File Upload option. The upper limit is currently 1000 sequences containing 3000 nucleotides per sequence.

The sequence in your question is a protein, so it consists of amino acids, not nucleotides. So, if you instead give a nucleotide sequence, for example this one:

>KX692212.1 HIV-1 isolate CH47 from Switzerland pol protein (pol) and gag protein (gag) genes, partial cds
CCTCAGATCACTCTTTGGCAACGACCCCTCGTCACAATAAAGATAGGGGGGCAACTAAAGGAAGCTCTAT
TAGATACAGGAGCAGATGATACAGTATTAGAAGAAATGAGTTTGCCAGGAAGATGGAAACCAAAAATGAT
AGGGGGAATTGGAGGTTTTATCAAAGTAAGACAGTATGATCAGATACTCATAGAAATCTGTGGACATAAA
GCTATAGGTACAGTATTAGTAGGACCTACACCTGTCAACATAATTGGAAGAAATCTGTTGACTCAGATTG
GTTGCACTTTAAATTTT

The program will work as expected:

For future reference, I found the KX692212 sequence above by taking your protein sequence and running a tBLASTn analysis against HIV sequences in NCBI's nr database. That was the first hit.


Finding simple sequence from reads with significant overlap I wrote a "script" to pull out reads from a huge fastq file in an iterative manner, by finding homology to the previous sequence.
It should be relatively easy to overlap them and assemble the gene, however traditional assemblers like SPAdes fail due to obviously low coverage.
I am hesitant to just feed it more faked reads, when I know I have all the data necessary to just overlap them into a single sequence.
The reads are set up like this:
GATCAAACATC
     AACATCAGTTAG
         TCAGTTAG
              TAGAGGATAGC

What kind of assembler is out there where I can feed it this kind of terribly-formatted data?

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


What does it mean to over sample? One of my studies (A) involve sequencing the microbiome. After selecting the variable region and the primers for the targeted region of the 16S sequence. The samples were sent to a platform.

For another study (B) we used the same strategy (same primers, region, machine...) but another platform. However we obtained from this platform a lot more sequences:
In study A the mean sequences per sample are 10k, while in the study B the mean sequences per sample are 80k

I consulted with others and it seems that the study A we under-sampled and in study B we over-sampled.

I can understand that by under-sampling we didn't obtain enough reads to be sure about the microbiome present in the sample. But what does it mean to over-sampling?

That expert I consulted suggested normalize the data to fewer sequences, discarding information we already have. I can understand that the sequencing has errors and can create chimeras but the resulting consensus sequences of the OTUs have found a match with databases of microorganisms. What or how should I take into account that I over-sampled?

As you sequence a library you are sampling from the pool of all possible oligonucleotides found in it. This is a stochastic process, meaning that you are more likely to sequence a piece of DNA that is more prevalent.

As you increase the number of reads that you sequence for a library, you decrease the chances of seeing a new read that you have not yet seen before. If you plot the number of unique reads per total reads sequenced, it will look something like the cumulative density function of the exponential distribution, where you have diminishing returns on unique sequences as you continue to sequence.

If you sequence a library and quickly stop seeing unique reads, we say that this library has low complexity. Therefore, you have "oversampled" it if you sequence too deeply. I am guessing that your metagenomics library actually has a lot fewer species/unique sequences than they anticipated, and that you sequenced deeply anticipating more.

Edit: In response to your comments, the bigger problem in 16S analyses is not sequencing too deeply, it is over-amplifying your samples such that they no longer represent what is actually in your gut microbiome samples. In this case if you used read normalization you would not accurately represent what was in the library that you sequenced since the amplified molecules (16S) are so small that you only have len(16S molecule)-(read len) + 1 possible reads in each direction. You might consider removing PCR duplicates with GATK but this will still have the same effect of not representing what was in the sequencing pool. It is very possible that you have many PCR duplicates after 16S amplification. I generally don't trust amplicon sequencing for quantitive analyses, and don't really trust for low coverage "is-this-species-there-or-not" either, especially when dealing with the 'dark matter' of uncharacterized species. You might be well off just asking "is this species present in high coverage at this time point?"


downloaded software in Conda environment does not show in ls command I installed gzip package in Conda environment, then  to make sure the software was there, I used "ls" command, it was not in the list. Tried to install again, system says all the packages have already installed, Why? How do I check the list of software installed in the environment? Thanks!

Just type
conda list gzip

OR
conda list | grep gzip

OR if its installed in a given environment
conda activate myenv

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


What programs account for structural alignment of different parts of distant homologs which have significant structural differences? If there is a need to perform structural alignment of different parts of distant homologs, which program one should use?
Since distant homologs often have significant structural changes, meaning the assumption of a similar physico-chemical environment for aligned residues would be violated.
Any program/tools or suggested reference will be a great help.
Thanks.

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


What samples can be used for a Mutect2 Panel of Normals I'm working on calling somatic variants from solid tissue tumors. I have one normal and one tumor sample for each patient. I plan to use Mutect2 to call somatic variants after preprocessing the data with the GATK Best Practices pipeline. Mutect2 recommends the use of a Panel of Normals (PON), but I can't figure out what samples should be used to built the PON. The guide says

  Additionally, the samples should come from subjects that were young and healthy to minimize the chance of using as normal a sample from someone who has an undiagnosed tumor. Normals are typically derived from blood samples.

My normal samples are from solid tissue, not blood.
Can I use the normal samples from each patient to build the panel of normals? Or would this cause problems because I'm using reusing the same normals for both the PON and as the normal sample for somatic variant calling? My first thought is that most somatic variants in individual normal samples should be filtered out by pooling all the normals together to build the PON but I wasn't sure if there are other reasons to avoid using non-blood based normal samples to build the PON

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

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


Removing Batch Effect in Heatmaps after Differential Gene Expression Analysis I'm working on a dataset in which the first replicate of each group is one batch and the second replicate is in a second batch. After checking the PCA plot and seeing the batch effect in PC1, I used removeBatchEffect function from limma to subtract the batch effect from my count data. Then, using PCA gives me a plot that doesn't seem to have any apparent batch effect left! However, it is recommended not to use batch effect correction for differential gene analysis but use the batch variable along with the group variable in constructing the model.matrix. So, I did that, ran limma/voom on the normalized counts, and extracted differentially expressed genes. However, when I'm trying to make a heatmap from the DEGs, I still see that the samples from different batches are clustered separately, instead of seeing the clustering of replicates of the same sample. So, my question is whether I should use removeBatchEffect on the count data from DEGs and then use the transformed dataset for heatmaps or there is another way to fix this?

It is indeed true that for the DE analysis one should include batch into the formula to avoid changing the original counts. Still, for everything else such as plotting heatmaps use of removeBatchEffects is perfectly fine and (at least to me) a standard and well-accepted procedure. It essentially does not matter what you use to correct for the batch effect for the counts you use downstream. Results will probably be similar. Combat-Seq from the sva package is a recent adaptation of ComBat specifically for RNA-seq which (from what I understand) better deals with the integer-count nature of the data. This operates on raw counts and avoids the infamous negative values that happen at times with both limma and Combat. After applying ComBat-Seq on your raw counts you can normalize them as usual with edgeR (or any tool you like) and then make the heatmaps. See https://github.com/zhangyuqing/ComBat-seq. The DE results still should come from the normal DE pipeline with batch as covariate as discussed above.


How to extract / convert gff3 CDS sequences to multifasta I would like to extract all the CDS from a batch of genomes. I have found a perl script from BioStars but this does not seem to work for me. I would preferably like to have a script/ method which will use the locus tag as a header.

Example

My gff files are from prokka output

ATCC0000    Prodigal:2.6    CDS 243 434 .   +   0   ID=FKKLIMLP_00001;Parent=FKKLIMLP_00001_gene;inference=ab initio prediction:Prodigal:2.6;locus_tag=FKKLIMLP_00001;product=hypothetical protein
ATCC0000    prokka  gene    243 434 .   +   .   ID=FKKLIMLP_00001_gene;locus_tag=FKKLIMLP_00001
ATCC0000    Prodigal:2.6    CDS 1727    2131    .   -   0   ID=FKKLIMLP_00002;Parent=FKKLIMLP_00002_gene;inference=ab initio prediction:Prodigal:2.6;locus_tag=FKKLIMLP_00002;product=hypothetical protein

I like bedtools getfasta.  My typical option set is bedtools getfasta -fi <reference> -bed <gff_file> -name -s.  Be aware of the -s to make sure you are pulling the correct strand.  I like bedtools because it is a versatile tool overall for handling bed, gff and vcf file manipulations.

# bedtools getfasta

Tool:    bedtools getfasta (aka fastaFromBed)
Version: v2.26.0
Summary: Extract DNA sequences from a fasta file based on feature coordinates.

Usage:   bedtools getfasta [OPTIONS] -fi <fasta> -bed <bed/gff/vcf>

Options:
    -fi Input FASTA file
    -bed    BED/GFF/VCF file of ranges to extract from -fi
    -name   Use the name field for the FASTA header
    -split  given BED12 fmt., extract and concatenate the sequencesfrom the BED "blocks" (e.g., exons)
    -tab    Write output in TAB delimited format.
        - Default is FASTA format.

    -s  Force strandedness. If the feature occupies the antisense,
        strand, the sequence will be reverse complemented.
        - By default, strand information is ignored.

    -fullHeader Use full fasta header.
        - By default, only the word before the first space or tab is used.


Loop Pymol command get_area I want to select multiple residues from different PDB files loaded on the same pymol session.
select resi LIG
get_area sele

How do I loop the above command to enable the command below:
select ligands,resn LIG

from https://www.blopig.com/blog/2019/08/how-to-iterate-in-pymol/
Example 1: Print the names of all the selected objects
You’ve selected several objects in PyMOL, and you want to see their names. By default, a new selection object called “(sele)” is created. We can create a Python variable called seleobjs, say, and then loop over them:
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: print obj

Example 2: Disable all the selected objects
If you still have your selection, “(sele)”, you can disable them (effectively turn them off on the right-hand sidebar):
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.disable(obj)

Example 3: Delete all the selected objects
You’re done with the objects in your selection, “(sele)”: you might be able to delete them using “delete ”, but if you can’t write a wildcard that covers all the names of the objects, you can iterate over a selection and delete them:
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.delete(obj)

from https://pymolwiki.org/index.php/Get_Area
get_area calculates the surface area in square Angstroms of the selection given. Note that the accessibility is assessed in the context of the object(s) that the selection is part of. So, to get the surface areas of e.g. a component of a complex, you should make a new object containing a copy of just that component and calculate its area.
please let us know if it worked


Loss of predictive power of polygenic risk score when dataset contains missing variants I am trying to calculate polygenic risk scores (PRS) scores for a new dataset. This dataset does not have all the variants that the PRS score needs. The PRS score I am interested in has 40 variants, but in my dataset only 20 of those variants are present. What is the correct way of assessing the decreased predictive power of the PRS score due to the missing variants?
I assume it is not as easy as summing the effective weights of those variants I have and comparing to the full sum (all variants needed for the PRS).
Thanks a lot and sorry for newbie question!

Which software package around you using, PRSice? Generally missing values comprising half the total observations are extremely difficult to calculate using regression analysis (or many statistical analyses for that matter.
One way forward is to include (an) example(s) of all 40 variants into your data set and perform the analysis with and without the added data, then compare this against a data set where all 40 variants are observed and your samples are excluded. You might find the missing values have low probability of PRS and this would then explain your results.
Missing analysis in regression is a specialist area of statistics and is not for the faint hearted, so a work-around is the best way forward.


How to find correlation between two specific genes in same dataset? I would like to plot the correlation between two specific genes in my data. I have a matrix with genes in rows and samples in columns, with read counts data. If I want to look at the correlation between 2 genes

1) do I need to convert the counts to some other units?

2) Is normalization needed?

I currently have counts returned by featureCounts:

          Sample1  Sample2  Sample3  Sample4  Sample5

ESAM         803    338       528       841      712
ESCO1        1594   1401      847       392      611
ESCO2         274   778       288       77       204
ESD          6192   3502      2013      1573     1808
ESF1         1356   1497      1011      372      701
ESM1         398    58        114       582      135

I tried this way:

I transposed the above data making genes as columns and sample names as rows. I removed the rows. Lets say the genes as columns is in dataframe "df".

df.mat <- cor(df)
colnames(df.mat) <- NULL
rownames(df.mat) <- NULL
heatmap.2(df.mat, trace="none")

I have an error:

Error in .External.graphics(C_layout, num.rows, num.cols, mat, as.integer(num.figures),  :
  invalid graphics state

I want to look at the correlation between ESCO1 and ESCO2. Can anyone tell me how to do that and plot with a small example.

I suppose in an idealized world you'd use something like salmon to get TPMs, which can account for differences in effective gene length by sample. However that's quite likely to be overkill.
Not necessary, you're using a with-sample ratio, which wouldn't be changed by normalization.

Here is a small example in R, assuming your matrix or data frame is named m:

first = which(row.names(m) == "ESCO1")
second = which(row.names(m) == "ESCO2")
plot(m[first,], m[second,])

or cor(m[first, ], m[second,]). You'll want a rank-based correlation.


Where can I find the DNA sequence of the C. elegans rRNA cluster? Where can I find the DNA sequence of the C. elegans rRNA cluster?

Ideally, I would be able to do it using the command-line.

I've tried using NCBI but the results are unsatisfactory. I am interested in the DNA sequence that is in the current reference genome of WormBase.

I have tried using the Bioconductor biomaRt package and create a table containing the description but when I grep 'rRNA', I only get proteins doing things with rRNA. I used the following code:

ensembl <- biomaRt::useMart('ensembl', dataset = 'celegans_gene_ensembl')
wormGenes <- biomaRt::getBM(attributes = c('external_gene_name', 'chromosome_name', 'start_position', 'end_position', 'gene_biotype', 'description'), mart = ensembl)
wormGenes[grep('rRNA', wormGenes$description), 6]

[1] "rRNA 2'-O-methyltransferase fibrillarin  [Source:UniProtKB/Swiss-Prot;Acc:Q22053]"
[2] "Pre-rRNA-processing protein pro-1  [Source:UniProtKB/Swiss-Prot;Acc:Q22006]"
[3] "Pre-rRNA-processing protein TSR1 homolog  [Source:UniProtKB/Swiss-Prot;Acc:Q19329]"
[4] "rRNA methyltransferase 2, mitochondrial  [Source:UniProtKB/Swiss-Prot;Acc:O62251]"
[5] "Probable rRNA-processing protein EBP2 homolog  [Source:UniProtKB/Swiss-Prot;Acc:Q09958]"
[6] "Putative rRNA methyltransferase  [Source:UniProtKB/TrEMBL;Acc:Q5CCJ2]"

The info is actually in the gene_biotype column of my wormGenes data frame. I could select the rRNA genes using dplyr::filter(wormGenes, gene_biotype == 'rRNA'). Unfortunately, no description is available for the genes, so I will have to get it using the gene IDs.


How do I create a for loop to filter through different FDR values? This is probably a really quick/simple fix but I am a Noob to R programming/programming in general. I am trying to create a for loop that runs the test multiple times with different FDR values (5,10,20). How do I make it run 3 times with different values and then produce 3 different datasets with names like... atac.glmtop.males.opening.fdr5, atac.glmtop.males.opening.fdr10, atac.glmtop.males.opening.fdr20? Also is there a way to make the write.table function work so that it exports files with the fdr5, fdr10, fdr20 in the name?

fdr.nums = c("5","10","20")
fdr.nums.vector = c(5,10,20)

for (i in 1:length(fdr.list)){
  atac.glmtop.males.opening.fdr[i] <- topTags(
    atac.glmfit.males,
    n = nrow(atac.dge$counts),
    sort.by = "none",
    adjust.method = "fdr"
  )$table %>%
  tibble::rownames_to_column("peakco") %>%
  filter(logFC < 0) %>%
  filter(FDR < (fdr.vector[i]/100)) %>%
  merge(.,annotated.atac.pbmc.narrow, by="peakco" ) %>%
  tibble::column_to_rownames("peakco")

  atac.glmtop.males.opening.fdr[i]

  write.table(
    atac.glmtop.males.fdr[i].opening[,c("chr", "start","end")],
    file = paste0("~/Downloads/sexM.age.groupHO.fdr",i, ".opening.bed"),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
  write.table(
    atac.glmtop.males.fdr[i].opening["GeneName"],
    file = "~/Downloads/sexM.age.groupHO.fdr"+ i + ".opening.GeneNames.txt",
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
}

For loops in most languages can be done in a variety of ways. What you have currently coded, loops over values 1 to 3 by using

for (i in 1:length(fdr.list)){
...
}

However, you probably want to use a different form and loop over each value in that vector using something more like:

for (fdr in fdr.nums.vector){
...
}

This way you have convenient access to the FDR value in each loop via the fdr variable. Which you can then use in constructing your output filenames. Thus your code would become something like:

fdr.nums.vector = c(5,10,20)

for (fdr in fdr.nums.vector){
  atac.glmtop.males.opening.fdr[fdr] <- topTags(
    atac.glmfit.males,
    n = nrow(atac.dge$counts),
    sort.by = "none",
    adjust.method = "fdr"
  )$table %>%
  tibble::rownames_to_column("peakco") %>%
  filter(logFC < 0) %>%
  filter(FDR < (fdr/100)) %>%
  merge(.,annotated.atac.pbmc.narrow, by="peakco" ) %>%
  tibble::column_to_rownames("peakco")

  atac.glmtop.males.opening.fdr[fdr]

  write.table(
    atac.glmtop.males.fdr[fdr].opening[,c("chr", "start","end")],
    file = paste0("~/Downloads/sexM.age.groupHO.fdr",fdr, ".opening.bed"),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
  write.table(
    atac.glmtop.males.fdr[fdr].opening["GeneName"],
    file = "~/Downloads/sexM.age.groupHO.fdr"+ fdr + ".opening.GeneNames.txt",
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
}

Advancing From for Loops to apply

As a beginner, for loops seem quite logical and easy to understand. However, R has a more powerful way to achieve the same thing through what is called "vectorisation" using the apply function or one of it's relatives lapply, sapply etc. Once you have understood and mastered for loops, you should take a look at apply, especially if you ever need to do a for loop many times or you find yourself doing nested for loops. R is notoriously slow with for loops and apply will likely be many times faster (potentially orders of magnitude). Here are some links for when you are ready to advance from for loops to vectorisation.

https://www.r-bloggers.com/how-to-write-the-first-for-loop-in-r/
https://www.r-bloggers.com/using-apply-sapply-lapply-in-r/


Biological replicates on Chip-seq Transcription factor data I have 2 biological replicates for chip-seq transcription factor data and I have applied macs2 peak calling for each replicate separately.

How can i make use of the biological replicates to extract the most out the data?

Thanks!

You could first look at the degree of correlation between the two replicates - what proportion of peaks are shared between the two, versus peaks found in one sample only?  This will give an idea of how repeatable the analysis is, and how many peaks are variable between samples or due to artefacts of the method.

When it comes to interpreting the biology behind the binding patterns, you are probably going to want to only interpret peaks which you find in both replicates.


Folded Protein Chunk Dimensional Classification? Are there known dimensional measurements for the classification of folded proteins given a starting chunk/domain as defined by something like the clustering functionality of MSM Builder?
Examples of what I would be looking for would be dimensions such as:

pH

Salinity

Dissolved Sugars

Reducing Agent Concentration

Ambient Electronegativity

With corresponding scales and precisions such as:

logorithmic vs linear

step size

distribution (e.g. dissolved/flat vs fractal/inclusions)

The above parameters would be included in folding software after a known starting conformation provided by MSM Builder or similar chunking software to derive a GUID for proteins without (or with extremely minimal and tuneable) collisions?
For instance, if hundreds of GUIDs on average end up equating to the same folded protein that's fine, but if one GUID could map to two or more thermodynamically stable folded proteins at greater than some incredibly small (e.g. 0.0000001%) then it wouldn't fit what I'm looking for.  I know this is inherently NP-hard, but the complete ID would be something like:
{
    sequence,
    chunkId,
    pH,
    salinity,
    dissolvedSugar,
    reducingAgentConcentration
}

So I would expect the chunkId and dimensional parameters to play a significant role in reducing collisions for different starting parameters.
It's worth noting I'm not seeking some magic bullet where you type an ID and get a fully folded protein, just a reliable system of GUIDs with enough starting information embedded in the ID to start the folding process.

CATH
The CATH database classifies protein by fold: https://www.cathdb.info/
So the value from that is probably the most useful for you.
Crystallographic conditions

pH
Salinity
Dissolved Sugars
Reducing Agent Concentration

Your crystallising conditions do not mean much. They are a solution that is close to precipitating your protein but slowly enough for it to crystallise. A protein can crystallise under multiple different conditions and a crystal may be chosen simply because it is a lower index well in the Hampton screening conditions...
These mean very very little. And unless you are planning on doing ML to determine the best screening condition these should be avoided.
Data from the PDB
You have probably noticed that the only stats PDB/PDBe give are crystallographic, like space group, unit cell dimensions, resolution, completeness in shells, Rfree, Rmerge etc. So nothing there is of much use... except for bound ligand. Often this is the cofactor, but often it's a cofactor analog or a product (as opposed to the substrate). But that information can come from Uniprot etc.
pH
A solution has a pH, a protein has a pI. At a pH equal to the pI it will be neutral —and less insoluble. ProtParam can calculate the pI from a linear sequence. However, surface charge totally shows different spots. APBS can map the solution to the Poisson-Boltzmann equation to your protein. There is a PyMOL wrapper, a Pyrosetta wrapper and more. But you still would need to some stats and 3D operations to convert it into a small vector of values —rototranslations and projects are not for the faint of heart or in a rush.
Melting temperature
BRENDA has lots of enzymatic data on enzymes including melting temperature or kinetics at different temperatures (which follow the Eyring equation (TST), until MMRT kicks in) which could be useful.
For a given protein you can calculate the relative folding ∆∆G (eg. with Rosetta Score), which is a proxy to melting temperature in a buffer it is happy in.
Despite paining me to say... some fanciful Poisson-Boltzmann map derived vector and ∆∆G are probably all you could get from analysing structures. Whereas the rest are from CATH classification and other databases that are not structure based.


Why am I getting empty expression data from GEO? I am trying to analyze the scRNAseq data from this study. In their Method section they write:

The accession number for the RNA and DNaseq data reported in this
paper is GEO:
GSE116237.

When I go ahead and try to pull this data in R using GEOquery, I am faced with a data frame with no features (i.e. no genes). The following code
library("GEOquery")
gse <- getGEO("GSE116237")
eset1 <- gse[[1]]
eset2 <- gse[[2]]

print(dim(eset1))
print(dim(eset2))

yields
Features  Samples
       0      184
Features  Samples
       0      681

I also didn't manage to download the data manually; I get an error message when trying to unzip the .zip files or get redirected to SRA Run Selector from which I also didn't manage to figure out how to download stuff.
Am I missing something or are these files simply corrupted? The paper is published in Cell, so I would assume this data is, to some extend, legit...

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


Getting a "system is computationally singular" error in sleuth I am analysing 142 samples belonging to 6 batches. Additionally, those samples belong to 72 strains, which means that for most of the strains there are two samples.

I could fit simple models (for strain and batches for instance), but when I get to the "full" model (~batch+strain), I get the following error:

so <- sleuth_fit(so, ~strain+batch, 'full')
Error in solve.default(t(X) %*% X) :
  system is computationally singular: reciprocal condition number = 5.2412e-19

I should point out that of the 72 strains, only 15 have samples in distinct batches. This means that most strains (57) have both samples in the same batch.

Is the error due to an unknown bug or rather to the experimental design? Does it mean that the information on batches cannot be used?

Thanks

EDIT I've posted the experimental design in a gist

batch   strain  replica
batch_1     strain_41   1
batch_4     strain_41   2
batch_1     strain_28   1
batch_4     strain_28   2
batch_1     strain_26   1
[...]

This happen when the variables (strain +batch) create a design matrix like this:

batch strain
1 1 #
1 1 #
1 2
2 2
3 3
4 3
...
16 72

Which means that some of the covariates are not linearly independent (ie batch 1 and strain 1), all the strain 1 is in batch 1.

You can correct for batch effects, but not in this design of the linear model (if you want to take into account the strain). You could do one batch more with those 15 strains that are in a single batch (if they are in different batch between them) that way you would get an independent design.

Many strains in a single batch are in the same batch. You need to increase the number of samples (recommended anyway due to the low number of samples per strain) to avoid this problem.

There are a lot of related question in Bioconductor support forum, from where I expanded an answer.


Merging two .txt files I have two .txt files and let's call it 'File A' and 'File B.Here are examples of what file A and file B looks like:

FILE A

#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno A
2   2000    2000    T   G   Pheno F, Pheno G
3   3000    3000    C   G   Pheno L

FILE B

#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno B
2   2000    2000    T   G   Pheno H, Pheno I
3   3000    3000    C   G   Pheno M
3   3000    3000    C   G   Pheno N

Here's how I want the output to look like:

#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno A, Pheno B
2   2000    2000    T   G   Pheno F, Pheno G, Pheno H, Pheno I
3   3000    3000    C   G   Pheno L, Pheno M, Pheno N

Both are tab-separated .txt files and I am trying to merge them together for SNP annotation and would love to hear suggestions on how to proceed this.

Whereas wouldn't be as fast as awk, here is a rather dirty R solution:

library(dplyr)
library(tidyr)
library(magrittr)

# File 1

# the "fill" parameter is used since not all the rows have the same number of cols
file_1 <- read.table("file_1.txt", skip = 1, fill = TRUE)

# paste() function with argument collapse is applied to the "phenotype"
# columns of each row (these cols start at position 6 and goes until
# the number of cols: dim(x)[2] gives the number of cols of a table
phenotypes <- apply(file_1[, 6:dim(file_1)[2]], 1, function(x) paste(x, collapse = ""))

# Cols referring to the genomic position, ref and alteration (cols 1:5)
# are combined with the newly created "phenotypes" column
file_1 <- data.frame(file_1[, 1:5], phenotypes)

# colnames are added, could be avoided if dealt with at file reading step
names(file_1) <- c("Chr", "Start", "End", "Ref", "Alt", "Phenotypes")

# File 2

file_2 <- read.table("file_2.txt", skip = 1, fill = TRUE)

phenotypes <-apply(file_2[, 6:dim(file_2)[2]], 1, function(x) paste(x, collapse = ""))

file_2 <- data.frame(file_2[, 1:5], phenotypes)

names(file_2) <- c("Chr", "Start", "End", "Ref", "Alt", "Phenotypes")

# The two tables generated as above are merged using the appropriate
# columns as "anchors" in the "by" argument of full_join()
# The resulting Phenotypes.x and Phenotypes.y cols are concatenated
merged_tables <- full_join(file_1,
                           file_2,
                           by = c("Chr", "Start", "End", "Ref", "Alt"))
merged_tables$Phenotypes <- paste(merged_tables$Phenotypes.x,
                                  merged_tables$Phenotypes.y,
                                  sep = ",")
merged_tables$Phenotypes.x <- NULL
merged_tables$Phenotypes.y <- NULL

# Duplicated rows (in terms of genomic pos, ref and alteration) in any
# of the files result in duplicated rows in the merged table and are
# "handled" by aggregating/summarizing the info at the phenotype cols
# by using paste()
# strsplit() is used to "split" the "phenotypes" on "," in order to
# remove repeating phenotypes resulting from the full_join() call

merged_tables %<>% group_by(Chr, Start, End, Ref, Alt) %>%
      summarize(Phenotypes = paste(Phenotypes, collapse = ",")) %>%
      ungroup() %>%
      mutate(Phenotypes = strsplit(Phenotypes, ","))

merged_tables$Phenotypes <- lapply(merged_tables$Phenotypes, function(x) paste(unique(unlist(x)), collapse = ","))
merged_tables$Phenotypes <- unlist(merged_tables$Phenotypes)

# A tibble: 3 x 6
    Chr Start   End Ref   Alt   Phenotypes
  <int> <int> <int> <fct> <fct> <chr>
1     1  1000  1000 A     T     PhenoA,PhenoB
2     2  2000  2000 T     G     PhenoF,PhenoG,PhenoH,PhenoI
3     3  3000  3000 C     G     PhenoL,PhenoM,PhenoN

# The resulting table (tibble in this case) can be saved as a file
# using write.table() or write.xlsx()

paste() in combination with its collapse parameter is used because of the unequal numbers of columns in the files. fread() of the data.table package could be helpful for reading large files.


How to remove batch effect from TCGA and GTEx data I want to play some differentially expressed genes for gliomas cancer. I collected cancer tissue profiles from the TCGA database and normal tissue profiles from the GTEx database. I straightforward use ComBat from R package "sva", and the result shows little differentially expressed genes (DEGs) were found. How can I remove batch effects correctly?
Thanks

This has been asked and answered here and Biostars. My take home message from many similar posts is "not to bother with correcting batch effects when performing differential expression analysis but to include 'batch' as a covariate". The latter accounts for the batch effects.
EDIT: This does not really answer OP's question, as @swbarnes2 and @ATpoint pointed out, they have a completely confounded design and the linear regression model would not accept such a model matrix.


Workflow of metabarcoding analysis For a diet analysis of an insect-eating animal, all species in a sample shall be identified. For this a sequencing of the metagenomic sample was done, where the COI/COX region was used as a barcode and amplified. This resulted in 20 million reads (after quality control).

While looking for an apropriate workflow, I am faced with these issues:

For many tools 20 million reads are to much to handle (BLAST for
instance)
For the majority of metagenomic tools the focus is on micro-organisms, not insects
The database to search against must contain COI/COX data, mainly from insects, however, the animal could also eat snails or frog eggs...
Since the reads are amplicons some analysis are not recommended

Leaving me with these questions:

What are the typical analysis steps in a workflow like this?
Which tools could I use for my purpose?
Where could I get a database of COI/COX data?

I'm working on something similar as a side project. I'll detail the steps that I'm using, but I'm curious if others can suggest a better way.

You don't mention if your data are paired-end, but in my case I have paired-end sequences that overlap in the middle. I used FLASh to merge the reads and seqtk to convert the reads to fasta format.
I used the cluster_fast algorithm in vsearch to cluster sequences at 100% identity, then filtered out rare sequence types (<10 sequences) to account for sequencing errors.
This resulted in a manageable number of unique sequences, which I searched against a custom database using vsearch, but which could also have been searched against the nr database using blast.


Tabix queries spanning chromosomes I can do something like this:

$ tabix ${url} ${range}

where ${range} spans a region on a single chromosome, e.g., chrA:X-Y.

Can a tabix query span multiple chromosomes, e.g., where ${range} looks like chrA:X-chrN:Y (or similar), which returns data between chrA:X and chrN:Y (using either a lexicographical or natural sort order between chrA and chrN)?

Note: I do not want to pass in a BED file of intervals. If it can't be done without passing in a separate file, that's fine. I only want to provide a one-liner to give to someone who can run tabix on their end, to do this type of query directly and remotely.

Looking at the htslib source code, it seems like the answer is no, a single region specifier cannot span multliple chromosomes.

Seems like the closest you can get is to enumerate all the chromosomes as separate regions on the command line. So to get everything between chrA:X and chrN:Y, you would have to write:

tabix inputfile.gz chrA:X chrB chrC ... chrN:1-Y

(If you omit the end coordinate like chrA:X, the region encompasses chrA from the specified coordinate to the end. If only specify the chromosome name without the coordinates, like chrB, the region encompasses the whole named chromosome, and for the last chromosome, you can specify up to the desired point by starting from the beginning: chrN:1-Y.)


SAM format: Does the BAM "Integer or numeric array" field no longer exist? Why? The SAM file format specification makes clear that optional fields must follow the format TAG:TYPE:VALUE

Type   Regexp matching                                       VALUE Description
A     [!-~]                                                  Printable character
i     [-+]?[0-9]+                                            Signed integer
f     [-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?                 Single-precision floating number
Z     [ !-~]*                                                Printable string, including space
H     ([0-9A-F][0-9A-F])*                                    Byte array in the Hex format10
B     [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+    Integer or numeric array

However, looking at the samtools source code, I noticed this:

https://github.com/samtools/samtools/blob/develop/sam_view.c#L920-L931

switch (aux_type) {
    case 'H':
    case 'Z':
        if (kputs(bam_aux2Z(s), linebuf) < 0) return false;
        break;
    case 'i': kputw(bam_aux2i(s), linebuf); break;
    case 'I': kputuw(bam_aux2i(s), linebuf); break;
    case 'A': kputc(bam_aux2A(s), linebuf); break;
    case 'f': kputd(bam_aux2f(s), linebuf); break;
    case 'B': kputs("*** Unhandled aux type ***", linebuf); return false;
    default:  kputs("*** Unknown aux type ***", linebuf); return false;
}

Questions:

(1) I have never seen an "Integer or numeric array" of case "B". What does such an tag look like? Is there an example?

(2) Why is this currently "*** Unhandled aux type ***"? Was this the case with previous versions of samtools?

An example would be ZC:B:i,271,271,1,0, which is from TMAP. The i part denotes the values held in the array, which are comma-separated.
That code is only used in samtools fastq, where that tag type isn't output into fastq files.


Ancestry of the coronavirus 2019-nCov, WuHan city, China In one of the answers to another question about the corona virus a link was given to this phylogenetic analysis of the virus.

Can somebody assist a non-bio type here? It seems to show that the current corona virus split from a virus in bats. And that the same ancestral virus is also ancestral to several other viruses, some in humans and some in bats.

Is my understanding basically correct?

Overview
The central focus of the tree is to highlight the key biological concern of the new coronavirus, 2019-nCov. The key concern is the genetic similarities to SARS epidemic, and relates to the SARS receptor.

SARS background SARS is endemic in bats (your BioRxiv tree partly shows that and this tree definitely shows it) and in the 2002 epidemic infected civet cats which then infected humans. More importantly it is then transmitted from one human to another. This is a consequence of the SARS receptor been able to exploit the cellular receptors in the respiratory tract of bats, civet cats and humans and use this tissue as a site of replication. The concern is exacerbated because the tree here shows SARS independently infected humans on two separate occasions, suggesting the cross-over has an underlying genetic basis, which is part of the "SARS lineage".

BioRxiv 2019-nCov tree Your tree shows the 2019-nCov has a recent common ancestor with SARS in comparison to the rest of the betacoronaviruses. This therefore provides some circumstantial evidence the receptor mechanism and ability to frequently crossover from bats and ultimately transmit human to human could be shared with SARS. At an evolutionary level (which is what ANY tree is supposed to reconstruct) it raises the question which common ancestor of SARS did the "SARS receptor" originate. I'd need to draw a diagram to better demonstrate this point, but I hope you get the idea.

The "influenza receptor" The analogy is with influenza virus and the entire epidemiology of influenza hinges around is sialic acid receptor of influenza and its ability to bind on to the cellular receptors in the upper/lower respiratory tract of birds, pigs and humans. The hypothesised mechanism is called the mixing vessel theory and is the classic epidemiological understanding of how new pathogenic influenza pandemics occur. If you replace "birds and pigs" with e.g. "bats, civet cats and humans" you get the idea why this could be scary. We don't the intemediate host of 2019-nCov ... but I speculate there must be one, unless eating Cov-infected bats is common in China.

Technical details of the BioRxiv tree The tree is a nicely diverse selection of the beta-coronaviruses. The authors have rooted the tree using the outgroups of delacoronaviruses and gammacoronaviruses, so it is a good robust selection of outgroups which can be used to correctly identify the direct of evolution of betacoronavirus divergence. In tree theory (its formal name is phylogenetics theory), extensive rooting is good and minimises artefacts.

One of the earliest members of the betacoronaviruses to diverge is MERS (Middle Eastern respiratory syndrome), which form a single "clade" (all viruses share a unique common ancestor) and these represent around 50% amino acid divergence from SARS/2019-nCov. The selection of MERS in the tree from both camels and humans looks good. The selection of other betacoronaviruses looks great, I wasn't aware of the "ruminant clade"  at all involving buffalo, cow etc .. infections and there was an associated human infection. There are loads of bat isolates of the betacoronaviruses throughout the tree, but we eventually arrive at the "SARS clade". The authors show that 2019-nCov is an outgroup to the SARS clade and shows a close relationship to one (BioRxiv tree) or two (this tree) bat isolates. Looking at the precise SARS clade (better shown in this tree) we can see loads of bat virus associating with SARS lineages. We therefore assume the reservoir to SARS and likely 2019-nCov is bat, moreover that the single ancestor to both viruses was a bat (it is called a parsimonious hypothesis).

The one thing the BioRXiv tree omits, for example with regards this tree, is the diversity of SARS and in particularly the two independent origins of SARS, which is a weakness of their analysis, particular if this fed into downstream analysis. It is not to say the authors were wrong, but it was uncool.

Word of caution my understanding is the divergence between SARS and 2010-nCov is around 15% and this is a quite a large amount of genetic divergence, even if both viruses share a most recent common-ancestor, form part of the same clade and receptor peptide motifs. Nevertheless there is sufficient amino acid divergence to generate notable differences in epidemiology, clinical symptom and transmission.

MERS, SARS and 2019-nCov mortality rates It is worth noting that within the betacoronaviruses the mortality rates between difference clinically important viruses is very different. MERS has a mortality rate of 40-60%, SARS is around 10% but 2019-nCov is <2.5%. The 2019-nCov mortality rate is still important given it has infected more and could potentially infect alot more people than SARS.


What are the applications of DNA or RNA pattern matching? I'm assuming we don't do pattern matching in DNA or RNA for the fun of it. So I'd like to know what are the applications of pattern matching or where does it fit in in larger applications?
I'm a computer scientist and when studying certain topics in CS we see a lot of references to applying algorithms to areas of biology or genetics one of which is applying text pattern matching algorithms to DNA to find sub sequences. Whenever I see those kinds of examples I think "Oh that's pretty cool" but then I don't see how it's applied in the real world for benefit and I can't find anything to suggest where these algorithms are applied in these fields.

One useful application of pattern matching is for DNA sequence motifs. As an example of a direct application, there is a class of proteins called transcription factors (TFs), which bind to DNA and regulate the expression of genes. There are many TFs. Each TF likes to bind to its particular conserved sequence of DNA, with some wiggle room for some parts of the sequence.
Algorithms for finding and scoring these patterns (and mismatches in patterns) is important for identifying regions where regulation occurs, which proteins are involved, and for answering biological questions related to gene regulation and basic biological functionality — and malfunction, in case of disease.


What tools can I use to look up my alleles, genotypes, or phenotypes from my sequenced DNA (WGS)? (I am not a Bioinformatics expert, please forgive and educate me if I've used any wrong terms or assumptions here)
I bought a "Whole Genome Sequence" kit, which gave me the following VCF files (this service did not include any analysis, just the raw data):
cnv.vcf
indel.vcf
snp.vcf
sv.vcf
My understanding is that these files show the differences of my DNA relative to a reference human genome (hg19 in my case)
What tools and resources can I use to find out which known alleles/variants of genes I have?
An example: Blood type
According to Wikipedia, there are 6 common alleles of the ABO gene, and which alleles I have determine my blood type.
Looking at my VCF results, I have these SNPs and INDELs within that gene:
<removed>

Given this, how would I figure out my blood type?
I understand that my data may not be 100% accurate, and that many phenotypes are the result of more than one gene, environmental factors, or are not well studied at all, but this should work for some simple things ie blood type, COMT, BHLHE41, CYP1A2, etc

If this is your own DNA and you are trying to analyse yourself you should talk to a genetic counselor. I personally think it's cool, but I deeply understand the potential consequences.
What you want is a variant annotation tool. There are many, but I normally write my own. Which is simply a case of pulling down a database and using pandas to match variants with annotation.
There are many databases. The biggest popular database is gnomad and the biggest disease database is clinvar. Others include omim, snpdb and 1k genomes. I'm not sure about blood sorry but I'm sure you can find one.


What i5 index should I use on the Illumina sample sheet for an unindexed p5 primer? I have an upcoming run on a HiSeq X and most of the libraries in the pool have both i5 and i7 indices. However, some of the libraries were made with the IS4 p5 oligo and it is unindexed.
The IS4 indexing oligo is from the following paper. The sequence is below.

Meyer, Matthias, and Martin Kircher. "Illumina sequencing library
preparation for highly multiplexed target capture and sequencing."
Cold Spring Harbor Protocols 2010.6 (2010): pdb-prot5448.

                                 miseq index
                              5' TCTTTCC 3'
                                 >>>>>>>
 5' AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTT 3'
    ^------------------^         ^--------------------^
    flowcell binding seq        adapter-binding seq

When we sequence these libraries dual-indexed on the MiSeq, the machine reads into the adapter and returns the index sequence TCTTTCC.
Question
What will the index sequence be when the HiSeq reads it?
The more fundamental question is: How exactly do MiSeq and HiSeq machines read index sequences?
side note
Here is what an index p5 and unindexed p5 molecule look like:

Update: Sept 2019 - After running into this problem again, an Illumina representative agreed with Devon Ryan and said that setting up a run for single indexing is the best option when there are unindexed p5 ends on one or more libraries.

Illumina machines read the two indices in separate index sequencing reactions where they have sequencing barcodes specific for indices. If you have no index, then you should either not bother with creating either index read or mask that during demultiplexing. If you're mixing indexed samples with an unindexed sample, you can't really rely on any particular sequence being returned by the machine (i.e., just use the "undetermined indices" fastq file). The unindexed adapters tend to not have the full index sequencing primer binding site on them, so reads without index generally don't return much signal. In practice, they're likely to adopt barcode sequences of neighboring clusters. We see that a lot when we spike in PhiX. It has no barcode but you tend to get a small amount assigned to each sample either due to index hopping or (more likely) bleed through of signal from neighboring clusters (or small amounts of those clusters walking between nanowells and then their minority signal dominating during index sequencing).


Biopython: resseq doesn't match pdb file I have a PDB file, and I need to extract its residue sequence numbers (resSeq's). Based on manual inspection of the first few lines of the PDB file (pasted below), I would think that resSeq's should begin with 22, 23. However, Biopython's PDB module suggests otherwise (output attached below as well). I wonder if it's a Biopython bug or if I have problems understanding the PDB format.

ATOM      1  N   GLY A  22      78.171  89.858  59.231  1.00 21.24           N
ATOM      2  CA  GLY A  22      79.174  88.827  58.999  1.00 20.87           C
ATOM      3  C   GLY A  22      80.438  89.415  58.391  1.00 21.89           C
ATOM      4  O   GLY A  22      80.362  90.202  57.440  1.00 23.18           O
ATOM      5  N   LEU A  23      81.588  89.069  58.972  1.00 21.51           N
ATOM      6  CA  LEU A  23      82.895  89.555  58.527  1.00 20.80           C
ATOM      7  C   LEU A  23      83.288  89.020  57.162  1.00 22.41           C
ATOM      8  O   LEU A  23      82.889  87.923  56.788  1.00 22.93           O
ATOM      9  CB  LEU A  23      83.973  89.232  59.560  1.00 20.97           C
ATOM     10  CG  LEU A  23      84.225  87.818  60.062  1.00 13.32           C
ATOM     11  CD1 LEU A  23      85.448  87.888  60.939  1.00 15.24           C
ATOM     12  CD2 LEU A  23      83.035  87.258  60.829  1.00 12.21           C

...
for i in chain:
    print i.get_full_id()

OUT:('pdb', 0, 'A', (' ', 2, ' '))
    ('pdb', 0, 'A', (' ', 3, ' '))
...

Let's check it. I copied your file fragment into a.pdb.

>>> from Bio import PDB
>>> for res in PDB.PDBParser().get_structure('pdb', 'a.pdb')[0].get_residues():
...   print(res.get_full_id())
...
('pdb', 0, 'A', (' ', 22, ' '))
('pdb', 0, 'A', (' ', 23, ' '))

So BioPython (ver. 1.66) gives 22 and 23 as expected.


Differential expression between two types of cells of specific marker gene expression I want to define the identity of some cells based on a gene marker (for example the colour yellow< EYFP gene). One group of cells has the EYFP gene the other doesn't and I want to find the differential gene expression between those groups that are found on the same cell sample. Do you have any tips for a code in Seurat? Thanks

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

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

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


Segment cell organelles with pixellib I've got some images of cell organelles and I really want to avoid labeling them by hand. The images all look something like the image below. Is there already an existing model specifially for cells? I tried it with existing models like mask_rcnn_coco.h5 or deeplabv3_xception65_ade20k.h5 but as they aren't built for cells the results aren't good.
My code is pretty basic right now:
from pixellib.instance import instance_segmentation

from PIL import Image

    def segmentImage():
        segment_image = instance_segmentation()

        segment_image.load_model("mask_rcnn_coco.h5")
        segment_image.segmentImage("cell.jpg", output_image_name="cell_seg.jpg", show_bboxes=True)

Cell example image:

Example segmentation result:

I am the creator of PixelLib. PixelLib supports instance segmentation of 80 classes of objects using pre-trained coco model. Cell organelles is not part of the classes supported.
The only solution is to train your own custom model. PixelLib supports custom training of a model using your own labelled dataset. Read this article I published on how to train a custom model with PixelLib.


GRCh38 centromeres mappability Is the mappability of the centromeres in the GRCh38 genome reference similar to each other?

As far as I can remember when GRCh38 came out, the sequence of the centromeres was determined by a combination of sequencing data and software prediction.

Given the way the sequence of the centromeres has been determined, should we expect Illumina 2x150bp (or shorter, 2x75bp) reads to map relatively equally to all centromere sequences?

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

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

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

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

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


Finding residues physically interacting in a PDB structure and sequence related to that structure I am going to describle the problem completely. First I should say that there are two kinds of interactions: transient and stable interactions. We are supposed to work on stable interaction like interactions between two monomers in a heterodimer. In a heterodimer there are two chains in which there are some residues between monomers which take part in the interaction of monomers and building the structure of heterodimer and also there are other residues interacting inside each monomer for building the structure of that monomer. There are some PDB IDs related to heterodimers in PDB. Now imagine that we are trying to show which residues are physically interacting between two monomers in a heterodimer and then we are trying two show the site of these physically interacting residues in the sequence related to that special structure. For that, we need to find the sequence related to that structure and we can download it from PDB. But if we want to make a multiple sequence alignment of that sequence, first we need to find the refseq sequence of that and then we must blast the sequence and find homologues of that sequence. Here there are some challenges: first, we must find residues interacting in the structure of heterodimer between two monomers. Second, we must find the sequences of monomer chains of the heterodimer and map the interacting residues in the structure to the sequence. Third, we know that when a structure is solved and its monomers are sequenced, may be they are not able to completely sequence the heterodimer and residue numbers are maybe different from the related refseq sequence. I would like to know how can we find the related refseq sequence (or sequences, in the case of heterodimer) and then how can we map the physically interacting residues in the structure of PDB ID to refseq sequence (finding the sites of these interacting residues in that refseq sequence) regarding three challenges mentioned above. Imagine that I have 100 PDB IDs of 100 nonredundant heterodimers and I would like to find and download the structures, sequence and finally refseq sequence. What should I do?

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

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

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

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

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

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

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

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

I hope this helps.


bedtool intersect multiple bed files I have 7 biological replicate of a normal sample after peak-calling I converted it into bed files and now I want to find what are the common region between those .

So i tried with bedtool intersect

my command is this
bedtools intersect -a SRR2920506.filtered.bed -b SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

it says " the intersect tool can detect overlaps between a single -a file and multiple -b files (instead of just one previously). One simply provides multiple -b files on the command line."

My doubt is the file that goes to "a" argument is not query which i have to detect overlaps against "b" files as these are all biological replicate. Am i doing the right way?

Secondly I tried with bedops but the number of overlaps which is reported is really low when i compare the result which i get from bedtools intersect function.

Any suggestion or help regarding how to intersect multiple would be really appreciated

If you want to find the intersection between ALL your bedfiles, you can try multiIntersectBed (available since bedtools 2.14.3).

It should work like this :

multiIntersectBed -i SRR2920506.filtered.bed SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

To sort your bedfiles by chromosome then by start :

sort -k 1,1 -k2,2n myfile.bed


What does these terms mean? I am working on cancer genomics failed to follow some terms even by googling

What is the difference of focal genome amplifications versus convergent amplification of given region of genome for instance fibroblast growth factors region?

Thank you for any clarification

I think these definitions still have a degree of subjectivity and I don't think they are mutually exclusive.

Convergent amplification basically means that "similar" cancers/neoplasms have "similar" amplification patterns, i.e. at around the same loci (see Myllykangas S et al.,2006 for an example).

Focal genome amplification is defined as "CNA [copy number alteration] of limited size, frequently enriched for cancer driver genes" (Krijgsmanthis et al., 2014).


ANY (technical) reason behind submitting sequences to GenBank versus ENA Sequence The DNA sequence sections of the three INSDC databases (i.e., DDBJ, ENA Sequence and GenBank) are synchronized periodically and strive to keep their stored data as ubiquitously accessible as possible. Except for idiosyncrasies in their data submission routes, there should be little, if any, reason for preferentially submitting sequence data to one database over another. Yet many researchers display very noticeable preferences!

When asked point blank, many colleagues emphasize that these preferences are not linked to the database-specific data submission or retrieval interfaces available. Yet they cannot provide a technical reason for their choice either.

Maybe researchers base their preference on geographic vicinity ("buy local"), but maybe there are genuine technical reasons for their preference (even if they are unaware of these reasons themselves).

Thus, my question is:

Can you think of any technical reason - however minute or seemingly insignificant - in preferentially submitting DNA sequences to NCBI's GenBank over ENA Sequence (or vice versa)? For example, there may be differences in data storage or accessibility that are relevant to you.

Put differently: Short of flipping a coin, why would you (as a bioinformatics-prone end user) select one database over the other for your data submission?

--

Edit 1: Many SE posters attempt to reassure me that I can submit my DNA sequences to either database and that "the data will be just fine there". Having submitted data to both GenBank and ENA Sequence for 10+ years myself, this is not what this question is about. Instead, this question is about carving out genuine technical differences from a user perspective.

Edit 2: Since no answer has been posted yet despite a bounty, I figure this is a hard question as specified. Thus, I will broaden the question and now also accept non-technical reasons for preferentially submitting DNA sequence data to either GenBank or ENA Sequence. (That being said, I still hope for someone listing a technical reason.)

Opinion: I prefer to use EBI/ENA because the turnaround time for "customer service" on editing metadata about a submission is faster than on NCBI. Inquiries to change metadata on NCBI can go unanswered for up to six months in my experience.

Answer to your question: No, there are no technical reasons that make it preferable to submit data through one website or the other. All entries end up being hosted by both databases anyway. The preference is up to the user to select whichever upload interface they prefer.


How can I install bioconductor-gviz and use it in jupyter notebook? I tried to install gviz in a conda environment, but that library seems to be incompatible to python and r. I tried to setup a clean environment using

conda create -n r -c conda-forge r-essentials jupyter

and then add the library with:

source activate r
conda install -c bioconda bioconductor-gviz

getting

  UnsatisfiableError: The following specifications were found to be in conflict:
    - atk

I remove atk, but now I get:

  UnsatisfiableError: The following specifications were found to be in conflict:
    - bioconductor-gviz
    - r-bindr

Does anyone manage to use gviz from within an jupyter notebook?

I also tried to install gviz from a running R-notebook with:

install.package('gviz')

  Warning in install.packages :
    package ‘gviz’ is not available (for R version 3.4.3)

Same when I try 'bioconductor-gviz'.

# Install the BiocManager package first
install.packages("BiocManager")
library("BiocManager")

# Then install the Gviz package
BiocManager::install("Gviz")


Can I run STAR without an annotation file? I wish to use Rascaf to scaffold a fragmented draft genome.
For this, I need to provide a BAM file of aligned RNA-seq reads and the draft genome.
So, I indexed the draft genome with STAR like this:
STAR --runMode genomeGenerate --genomeDir output/index --genomeFastaFiles draft.fa

Then, I tried to aligned the reads like this:
STAR --genomeDir output/index --readFilesIn reads.fastq.gz --readFilesCommand gunzip -c --outFileNamePrefix output/alignment --quantMode GeneCounts --outSAMunmapped None --outSAMtype BAM SortedByCoordinate --outSAMattrRGline ID:RG1 CN:yy \"DS: z z z\" SM:sample

The latter commands outputs the following error:
Transcriptome.cpp:51:Transcriptome: exiting because of *INPUT FILE* error: could not open input file exonGeTrInfo.tab
Solution: check that the file exists and you have read permission for this file
          SOLUTION: utilize --sjdbGTFfile /path/to/annotantions.gtf option at the genome generation step or mapping step

It suggests I should use a GTF file. But this is a draft genome of an individual for which no GTF is available. I could try to adapt a GTF from a close relative but for the scaffolding, it seems superfluous. Is there a way I can map with STAR without the GTF?
These commands were run on a GNU/Linux machine.

I don't think you can use the --quantMode GeneCounts option with no annotations. I think the error is trying to look for an exon file generated from the annotations to do the quantitation on. Remove that and I think it should work, as the manual specifically states that annotations are optional but highly recommended.


What does it mean for an amino acid to match with itself 99% of the time in a PAM1 matrix? In this video, the lecturer says:

  You see that phenylalanine is probably still phenylalanine because only one percent of amino acids have changed so roughly 99 percent of the time, phenylalanine is still phenylalanine. Now why is it not exactly 99 percent? Because it’s one percent for all 20 amino acids so some are greater than 99 percent and some are less than 99 percent. It averages out to 99 percent.

What on earth does he mean by this? I either do not have the background knowledge or am not thinking hard enough, or I need this to be explained in a way that is more understandable.

It means that the 99% of conservation of amino acids is just a mean value for all the amino acids. That there can be amino acids that are less conserved and some that are more conserved.

So from the 21 basic amino acids, arginine is more conserved than lysine in those alignments. But in general they are conserved the 99% of times.


What is the per-gene score in BUSCO output? The BUSCO _full_table.tsv output has scores per gene (an example bellow).

In the original publication, authors just mention a condition to satisfy a bitscore, but I could not find what the bitscore actually is. In the use guide they say

  If found to be complete, whether single-copy or duplicated, the BUSCO
  matches have scored within the expected range of scores ...

but again, I don't understand what the score actually represents. I am even not sure if greater the better or visa reverse. Does anybody know what the score is?

# BUSCO version is: 4.0.2
# The lineage dataset is: insecta_odb10 (Creation date: 2019-11-20, number of species: 75, number of BUSCOs: 1367)
# Busco id      Status  Sequence        Gene Start      Gene End        Score   Length
47at50557       Complete        scf1:20448-42284    20448   42284   4607.1  3554    https://www.orthodb.org/v10?query=47at50557     Immunoglobulin-like domain
547at50557      Complete        scf2:9780-23590    9780    23590   2743.8  1555    https://www.orthodb.org/v10?query=547at50557    protein unc-80 homolog
755at50557      Complete        scf3:9036-28083    9036    28083   3149.9  2206    https://www.orthodb.org/v10?query=755at50557    SRCR-like domain

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


How to remove sequence reads contains more than 2 X from multifasta file? I have 5000 protein sequences in one multifasta file. I found more reads have gaps as X in their reads. So, want to eliminate those reads completely (Whole protein seq) from the file. I am keeping filter criteria as if a read contains morethan 2 X ( continuesly or anywhere in the read) should be removed. Thanks in advance for your help.

The input sequence looks like this

>Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF
>Prot2
ANSTVKKKKLLLYYYSSSEERXXXXXXXXXXXFGHYFGHYFGHFYVHFGFYVHCEDYHF
>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

I want output Like this

>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

A Biopytnon solution that works for .gz files too :-). The script make a script filter_faa.py:

#!/usr/bin/env python3

import sys
from Bio import SeqIO
from mimetypes import guess_type
import gzip

filename = sys.argv[1]

encoding = guess_type(filename)[1]
_open = partial(gzip.open, mode='rt') if encoding == 'gzip' else open

with _open(filename) as ffile:
    for seq_record in SeqIO.parse(ffile, 'fasta'):
        if seq_record.seq.count('X') < 3:
            print(seq_record.format("fasta"))

Then usage is

python3 filter_faa.py my.faa > my_filtered.faa

The my_filtered.faa then will contain

>Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF

>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

The filtering condition is hard-coded. You could potentially parametrise it if you feel like later you might want to filter by a different number of unknown aa.