PCA plot shows big difference but not many differentially expressed genes are found I got a PCA plot of bulk RNA-seq experiment that looks the following way:

It was generated by the following code:

pcaData <- plotPCA(rld_sva, intgroup=c("Group"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  coord_fixed()

First sva correction was run to correct for batch effects and then rlog transformed values were plugged in to plotPCA function.

The first issue that catches the eye is that 100% of the variance is explained by just 2 dimensions. I am not sure what can one say about the data in this case. The second issue is that I get only around 20 differentially expressed genes by using DESeq2 analysis (log2FoldChange > 1, p_adj < 0.05). I know that we can not directly state that if there is a large difference on PCA there will be present a plenty of differentially expressed genes, but why is it not the case? Simple logic tells me that pca shows the difference between the samples in their gene expression, so I would expect seeing a plenty of differentially expressed genes.

You only have 4 samples total. I think it would be difficult to not have the PCA show big differences between the groups with so few points.

On the other hand, for differential expression, it is hard to get something to be statistically significant with only 2 replicates.


problem of "ordering in physical positions" phasing SNPs with Shapeit I'm trying to phase my data (whole genome resequencing SNPs) with Shapeit (from ped file) and I get this error message :

[33mERROR:[0m Wrong ordering in physical positions curr_pos=54 prev_pos=65

I found that it's due to the fact that my .map file is like this (4 columns):

Chromosome_number    SNP_ID   SNP_genetic_position    SNP_physical_position
Contig100131            .              0                        65
Contig100220            .              0                        54

And Shapeit only looks at SNP_physical_position but does not check if the SNPs are on the same contig or not and so it says it is not happy because position 54 is after 65!

How could I make it understand that the positions are not in the same contig and so there is no mistake? Must I run Shapeit per contig (I have a lot of them so I hope not)?

The manual suggests that each "chromosome" needs to have its own input file.

You have to split the dataset by chromosomes. ShapeIT can just phase one chromosome at a time. If you are working on GWA data, you can proceed as follows:

for chr in $(seq 1 22) ; do plink --file example --chr $chr --recode --out example_chr$chr ; done

You will obtain 44 files; example_chr1.ped, example_chr1.map, ..., example_chr22.ped, example_chr22.map that you can phase separately.


Find gene at position from gff or gbk file I have a VCF file with SNPs from a bacterial genome and want to find if the SNPs are located inside genes, is there some CLI-tool where you can pass a VCF file and a gff or gbk file and it returns the name of the genes?

Via BEDOPS:

$ gff2bed < annotations.gff > annotations.bed
$ vcf2bed < snps.vcf > snps.bed
$ bedmap --echo --echo-map-id-uniq snps.bed annotations.bed > answer.bed

This can be reduced to a one-liner if you're using bash:

$ bedmap --echo --echo-map-id-uniq <(vcf2bed < snps.vcf) <(gff2bed < annotations.gff) > answer.bed


Why do BLASTn and prokka not seem to be searching the whole fasta file? When I use blastn and prokka (I will detail exactly how I did so below) on a 2.8 million bp fasta file I get output start/end numbers that do not seem to cover the entire genome.

Starting with a .fna genome such as genome.fna I ...

1 - blastn

Searched the genome for 3135 different 28-mers using BLASTn.

makeblastdb -in genome.fna -dbtype nucl -parse_seqids -out ./output/genome

next command in python

blast_tsv = NcbiblastnCommandline(query=Q, db=DB, perc_identity=100, outfmt=6, out=(OUT))
    stdout, stderr = blast_tsv()

Q is this list of k-mers. DB is the database created

This outputs a list of search results but the important thing is that none of the start/end (columns 7 and 8) integers were greater than 100,000 yet the entire genome is 2.8 million base pairs long. I can provide the link to this file in the comments. Does this have to do with blastn stopping searching after finding one match? And if so how can I tell blastn to search for every match for each k-mer? (I'm open to using other blastn programs other than biopython)

2 - prokka

I used prokka to create gff files for the genome.

prokka-1.12/bin/prokka —setupdb
prokka-1.12/bin/prokka -kingdom Bacteria -rfam -outdir ./prokka_database/genome -force -prefix "${genome/.fna/}" genome.fna

Th gff file produced by this command only described genes with start and end positions less than 200k. I can provide the link to this file as well.

With a k-mer size of 28 it shouldn't be finding that many matches. And the prokka results are suspicious as well. Maybe you have multiple contigs (none larger than 100kb) in that file? What is the result of

grep ^'>' fasta_file | wc -l

? This would show how many contigs you have in the file.


Where can I find a list of SNP rs IDs that belong to the X chromosome? Where can I find a list of SNP rs IDs that belong to the X chromosome?

I'm trying to exclude the SNPs on X-chromosome in a data set, but I don't know where these SNPs come from.

I tried to use dbSNP to search, it returns 26933249 SNPs, but the webpage seems to refuse me to download the information.

You just searched for X[Chromosome], you didn't specify a species. Presumably, your dataset comes from a specific species, so you should limit your search to that species. For instance, for human, you can search dbSNP for X[Chromosome] human[organism]. Then, click on the send to link at the top left and download the data:

This is veeeeeery slow, however. I suggest that you instead download the data for chrX from dbSNP's FTP site (this assumes your species of interest is human): ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/chr_rpts/chr_X.txt.gz

Once you have this file, you can use it to filter your dataset:

Collect all rsIDs from the file:

zcat chr_X.txt.gz | awk '/^[0-9]/{print "rs"$1}' > chrX.rsIDs

Filter your dataset.

You haven't shown the format of your dataset, so if it is just a list of rsIDs like this:

$ cat dataset
rs77669712
rs199706910
rs2799940
rs12205639
rs10958205
rs200778574

You can remove any that are on chrX with:

$ grep -vxFf chrX.rsIDs dataset
rs77669712
rs199706910
rs12205639
rs10958205
rs200778574

That will print all rsIDs from the text file dataset that are not also in the text file chrX.rsIDs. Note the -x option for grep which tells it to only consider matches that cover the entire line. This ensures that rs1234 is not removed when you're looking for rs123.


COG Annotation - Dealing with genes assigned to two or more COG categories I am dealing with a list of genes which have been selected from a gene enrichment analysis. In order to see what kind of genes are overrepresented, I ran eggnog-mapper to do an orthology assignment of each gene against a bacterial database (provided in eggnog).

After running eggnog you get a .csv file containing for each gene which is the closest ortholog, and a list of GO, KEGG and COG terms associated.

From there, I would like to do a Fisher-test (also considering Bonferroni correction) to highlight COG categories overrepresented in our dataset compared to our reference genome.

Finally my question: How do you deal with genes assigned to multiple COG categories (same would apply for GO terms). You count +1 to each of the COG categories assigned? Do you just choose one category per gene?

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

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


python does not quit when the input file is too big I have a script that extracts the set of unique kmer pairs from a kmer dump file (& then using them downstream to guess ploidy).

The script seems to work correctly (it's reporting a meaningful set of kmer pairs), but since we have added a parallelization it does not quit when the job is done for really big data sets. The problem is that I never managed to reproduce the behavior with a small reproducible example and it takes around 400G of memory and a day of computation to get to the problem. My imagination is quite limited, I have no idea how should I debug such a problem.

I don't really get it; the script actually prints the last line

  *_families_2.tsv and *_coverages_2.tsv files saved.

, but then it runs further (& forever). I thought it could be related to longer time to free the memory as discussed here, but it does not seem like the case because it does not quit even after a really long time (I tried to just leave it for a day or so and when i kill it the memory is freed in like 20s).

Does anybody an idea where the problem could be? Or how to figure out?

The reported problem on GitHub.

You've hit the RAM bottleneck. Its very common in de novo assemblies and resolved using template assemblies. I don't know this program but suspect it is performing a pairwise analysis for every piece of data in a genomics data set.
The solution is to either:

find a bigger machine, e.g. a big multi-core machine and make sure the program distributes across every core.

parallelise the calculation (an excellent way to overcome RAM problems), or,

find an approximation to the calculation that is less computationally demanding, templates are a good way to circumvent this.


how to Know the longest ORF in Protein? I have a protein Seq. and I want to know

how to Know the longest ORF in this Protein?

and if there is any tutorial for this process ??

Proteins don't have ORFs, DNA sequences do.  What DNA sequences do you have?  If cDNA (i.e. transcript/mRNA) then an ORF finder will give multiple ORFs and the coordinates on the full-length input DNA which indicate it's length.


How to calculate various properties of a protein structure per atom using PDB2PQR and Rosetta tools (or any other tools)? Currently, I am implementing a descriptor of protein structure and would like to calculate properties per atoms like charges, hydrophobicity, hydrogen bond donors / acceptors, hydrophobicity, aromaticity etc.
I tried using PDB2PQR and succeeded with protein structures but it crashed while calculating nucleic acid bound proteins and choosing other force fields did not help.
Also, as far as I know, Rosetta tool can calculate hydrophobicity however, it assigns the property per residue but not per atom (correct me if I am wrong).
Any other suitable tool suggestions would be appreciated.

Per Atom
You are correct. Rosetta scorefunction does not store any per atom data. The scoring operates at the per residue level. Whereas each atom has its coordinates and properties in full atom mode, in Pyrosetta it is clear that an atom is a just part of a residue and every operation is applied at the residue level. It's a team effort: the functional group of the side chain works as one and in most cases two atoms have the same role (e.g. OE1 and OE2 in glutamate behave identically).
So, the real question is whether you really care about the atom level in a classical mechanical setting?
Formal charges
A wee parenthesis, we need to make a difference between formal charge and Marsilli-Gasteiger partial charge —one is an integer and the other is a float. Both require a stated protonation state in most applications —say aspartate vs. aspartic acid.
PDB2PQR formats to a PQR file for the ligand docking application autodock. It does not assign how a residue is protonated.
In the case of Rosetta, the partial charge depends on the 4 column of the topology (params file) of the amino acid, i.e. it is fixed for all amino acids —Rosetta does not implement a Drude particle to make something polarisable. I very strongly suspect this is the same in PDB2PQR.
Therefore, something like the Advance Poisson-Boltzmann Solver may give more interesting data in terms of partial charges, but also requires the uses to specify the protonation state.
A small parenthesis is histidine, which can adopt four forms. Two tautomers at neutral pH: single proton on either the epsilon (HIE in Amber) or delta nitrogen (HID), while it can be deprotonated (anion) or diprotonated (HID, cation). Rosetta switches between these two automatically.
Properties
You list a few properties of interest, which all seem to be based on the residue a priori of structure.
For example:

hydrophobicity, unless you have a polarisable model (Drude particle) this is structure independent. RSA, relative surface area, is a must better indication (per residue), which then can get classified as core or surface based on this. Per atom would require just a wee tinker.
aromaticity, this is also structure independent, but π-π stacking or π-cation interactions are noted at the residue level.
hydrogen bond acceptor/donor. Whether an atom is a donor or acceptor is actually defined in the atom type used in the topology file. Whether that atom fulfils that bonding is a different matter. If you really really wanted to get per atom values of this, in Rosetta you could exclude the backbone atoms you could mutate to alanine without minimising and see if it's maintained. So differentiating between backbone and sidechain is fine, but per atom is a different matter.

A structure-dependent property is how flexible the residue is, as measured as an RMS. This can be calculated per atom. This is akin to B-factors, which is how big is the "blur" of the atom. However, if you get the NZ atom on a lysine that is switching between two or more positions you get a partial occupancy. If this atom is basically hand-waving along an axis in an arc you get an anisotropic B-factor. In other words it does not vibrate in all 3 cartesian axes equally. So a B-factor is a terrible approximation. An RSMD of an ensemble is even worse.
In Rosetta, you can a quick vibration using Backrub, which is very quick at telling you how much motion does an atom experience.
Here be dragons

Warning: As mentioned, it takes two atoms to resulting in a force, hence why per atom scores are a bad idea.

However if for some exceptional reason this is required, say
figuring out which group on a ligand is causing a clash, in pyrosetta it is possible:
import pandas as pd
import pyrosetta
pyrosetta.init()
pose = pyrosetta.pose_from_sequence("ELVISISALIVE")
target_resi = 9
residue = pose.residue(target_resi)
scores = []
for i in range(1, residue.natoms() + 1):
    iname = residue.atom_name(i)
    for r in range(1, pose.total_residue() + 1):
        other = pose.residue(r)
        for o in range(1, other.natoms() + 1):
            oname = other.atom_name(o)
            score = pyrosetta.toolbox.atom_pair_energy.etable_atom_pair_energies(residue,
                                                                                i,
                                                                                other,
                                                                                o,
                                                                                scorefxn)
        scores.append(dict(zip(['residue_index_A', 'atom_name_A', 'residue_index_A', 'atom_name_B', 'lj_atr', 'lj_rep', 'fa_solv', 'fa_elec'],
                         [target_resi, residue.atom_name(i), r, other.atom_name(o), *score])))

df = pd.DataFrame(scores)
df.pivot_table(index='atom_name_A',
               aggfunc='sum')[['lj_atr','lj_rep','fa_elec','fa_solv']]


I have really skewed RNA-seq data, what's the best way to normalise it? Preferably in R! My data frame compares the RNA-seq reads from many genes in different tissues.
The reads look as follows.

I tried using log to make it better but still looks pretty skewed.

Are there any better transformations/functions for this in R?
This is part of an assignment I have, I have used in the past edge software and agree it will do a better job but in this work I'm not meant to use any of the pipelines.

Welcome to RNA-seq. This is not at all unusual, as RNA-seq raw data are not normally distributed. People typically approximate it via the Negative Binomial distribution. I suggest you check existing packages as Devon suggests in his comment for normalization and differential expression (or anything else). Typical choices are DESeq2, edgeR or limma, but there are others as well.
This guide might help you get started. Please also use google to explore existing threads both here at SE and over at biostars.org.
Most basic questions have already been discussed there and are informative to get started.


How can I locate duplicated regions in a sequence? I am facing an issue when trying to align short reads against a region in human chr5. The two Sensory Motor Neuron genes, (SMN1 and SMN2) are almost 100% identical and this causes the aligner to fail to align reads correctly since each read matches perfectly to two separate regions of the reference, leading to very low alignment scores.

If I BLAT the SMN1 gene's full genomic sequence against the human genome (hg19) on UCSC, I get (only showing the 1st two lines):

QUERY       SCORE START   END QSIZE IDENTITY  CHROM  STRAND  START       END   SPAN
--------------------------------------------------------------------------------------------------------------
NG_008691.1 35072     1 35072 35072   100.0%  chr5   +    70215768  70250839  35072
NG_008691.1 34279     1 35072 35072    99.8%  chr5   +    69340350  69375383  35034

And on the browser:

To get around this, I want to mask the SMN2 gene. However, the region of similarity extends beyond the annotated boundaries of the gene, so I would like to identify the full extent of this duplication.

So, given a large nucleotide sequence in fasta format, how can I compare it to itself and identify identical or near identical regions?

If you know approximately where are the duplicated regions, you can cut them out from the genome and align the two block using tools like MUMmer or Mauve.

I downloaded the sequence of ch5 you linked and extracted big regions around the genes (300bp blocks). The first block is from 70100000 - 70400000 (with SMN1 sequence at 115768 - 150839) and the second from 69200000 - 69500000 (with SMN2 sequence at 140350 - 175383).

I used progressive mauve to align the two (click on File -> align with ProgressiveMauve -> add Sequence ... -> Align).

I found that the duplicated region is almost 200k long and maps between 44690 - 225520 of the first block and 19760 - 200940 of the second block, which we still need translate to the original genome coordinates, but I believe I can leave this to you :-).

P.S. you can also somehow upload annotations to Mauve, but I was lazy, so I "hand drawn" the arrow where I knew the genes are.


load the phenotype data for ballgown I am trying to reproduce the work of this paper
[1], and I have run StringTie successfully, but after that I have to run Ballgown but could not understand this command:

>pheno_data = read.csv("geuvadis_phenodata.csv")

Actually, I dont know how to make the .csv file of my samples. What phenotype information should this csv contain?

I have checked the supplementary data of this research article but didn't find any .csv file.

How I can make the .csv file for my samples to load the phenotype data for ballgown?

  Load the phenotype data for the samples. An example file called
  geuvadis_phenodata.csv is included with the data files for this protocol
  (ChrX_data). In general, you will have to create this file yourself. It contains
  information about your RNA sequencing samples, formatted as illustrated in
  this csv (comma-separated values) file. Each sample should be described on
  one row of the file and each column should contain one variable. To read this
  file into R we use the command read.csv. In this file the values are
  separated by commas, but if the file were tabdelimited you could use the
  function read.table.

    pheno_data = read.csv("geuvadis_phenodata.csv")

Pertea et al, Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown, Nature Protocols volume 11, pages 1650–1667 (2016)

You can easily do it via Notepad++ or even Excel. Just don't forget to save the output file as a .csv file (comma-separated value).
The info inside should look like (my working example):

id,state
lar_control141,control
lar_control142,control
lar_control143,control
lar_probiotic142,probiotic
lar_probiotic143,probiotic
lar_probiotic144,probiotic
lar_probiotic145,probiotic
lar_probiotic146,probiotic

That should work


How to calculate overlapping genes between two genome annotation versions I have two annotations of the same genome generated with different annotation pipelines. I want to identify overlapping gene models.

An important feature of this genome is that there are many 'genes within genes', i.e. a genemodel in the intron of another genemodel. Therefore, I only want to count two genemodels as overlapping when their coding sequence exon annotations overlap.

Using something like bedtools intersect it is straightforward to calculate overlap between the gene-level annotations.

However: I am not sure how to select genes as overlapping when only their coding sequence exons (CDS features) overlap.

Short Answer:
In my opinion, my approach would be to pull out the CDS exons and run bedtools on those.

A Few More Details:
When you pull out the exons, make sure that you assign them all IDs if the don't already have them assigned and record which IDs "belong" to which genes. Now when you get exons that overlap, you know that they are coding and you can tie them back to which genes they originate from.


Blast results filters I performed a blastn search of NGS data against ssRNA database download from Internet, with a expected value 10-4. The size of NGS data reads is of 125 bp.
I have analyzed the blast results of the libraries. These libraries are from different habitats. The most difficult habitat, oakwood, produced a blast results with a maximum length of 51 bp. These results have a qcov of 100 (in most of the cases).
I would like to know if these low length sizes comes from degradation?
Should do  I filter these samples with a length = 125 bp with the aim to obtain results with better quality or with a qcov of 100 is enought to get reliable results?
Thanks you in advance.

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

!


How to reduce the occupied RAM when you are dealing with a very sparse matrix in a single-cell Experiment in R? I'm dealing with a very large and sparse dataset and the first issues I met occurred when I tried to use quickCluster that reported me this error:
                'cannot allocate vector of size 156.6 Mb'

So, given that I cannot wait to change the RAM of my computer and I can't afford to use a cluster, I want to rely on some other strategies like some package that would allow me to handle sparse matrices. I'm thinking about sparseM but given that I don't know well this package I'd like to know how to shrink the ram allocation for these kind of matrices. Any suggestion will be very appreciated!
Link to the dataset

Ah, looks like I can't even procrastinate on StackExchange anymore without seeing work-related stuff. Oh well.
Anyway, the other answers and comments are way off. scran has supported sparse matrices for years, ever since we switched over to the SingleCellExperiment class as our basic data structure. quickCluster does no coercion to dense format unless you tell it so explicitly, e.g., with use.ranks=TRUE (in which case you're asking for ranks, so there's little choice but to collapse to a dense matrix).
You don't provide an MWE or your session information, but this is how it rolls for me:
# Using the raw counts in the linked dataset. Despite being
# called a CSV, it's actually space delimited... typical.
library(scater)
mat <- readSparseCounts("GBM_raw_gene_counts.csv", sep=" ")

# Making an SCE just for fun. Not strictly necessary for
# this example, but you'll find it useful later.
sce <- SingleCellExperiment(list(counts=mat))

library(scran)
system.time(clust <- quickCluster(sce))
##   user  system elapsed
##  3.170   0.174   3.411

This is running on my laptop - 16 GB RAM but I'm definitely not using all of it. I only go full throttle when I'm working on some real data, e.g., the 300k HCA bone marrow dataset. Check out the book for more details.
Session info below, I don't know quite enough SO-fu to know collapse it.
R version 4.0.0 Patched (2020-04-27 r78316)
Platform: x86_64-apple-darwin17.7.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS:   /Users/luna/Software/R/R-4-0-branch/lib/libRblas.dylib
LAPACK: /Users/luna/Software/R/R-4-0-branch/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] scran_1.16.0                scater_1.16.1
 [3] ggplot2_3.3.2               SingleCellExperiment_1.10.1
 [5] SummarizedExperiment_1.18.1 DelayedArray_0.14.0
 [7] matrixStats_0.56.0          Biobase_2.48.0
 [9] GenomicRanges_1.40.0        GenomeInfoDb_1.24.2
[11] IRanges_2.22.2              S4Vectors_0.26.1
[13] BiocGenerics_0.34.0

loaded via a namespace (and not attached):
 [1] beeswarm_0.2.3            statmod_1.4.34
 [3] tidyselect_1.1.0          locfit_1.5-9.4
 [5] purrr_0.3.4               BiocSingular_1.4.0
 [7] lattice_0.20-41           colorspace_1.4-1
 [9] vctrs_0.3.1               generics_0.0.2
[11] viridisLite_0.3.0         rlang_0.4.6
[13] pillar_1.4.4              glue_1.4.1
[15] withr_2.2.0               BiocParallel_1.22.0
[17] dqrng_0.2.1               GenomeInfoDbData_1.2.3
[19] lifecycle_0.2.0           zlibbioc_1.34.0
[21] munsell_0.5.0             gtable_0.3.0
[23] rsvd_1.0.3                vipor_0.4.5
[25] irlba_2.3.3               BiocNeighbors_1.6.0
[27] Rcpp_1.0.4.6              edgeR_3.30.3
[29] scales_1.1.1              limma_3.44.3
[31] XVector_0.28.0            gridExtra_2.3
[33] dplyr_1.0.0               grid_4.0.0
[35] tools_4.0.0               bitops_1.0-6
[37] magrittr_1.5              RCurl_1.98-1.2
[39] tibble_3.0.1              crayon_1.3.4
[41] pkgconfig_2.0.3           ellipsis_0.3.1
[43] Matrix_1.2-18             DelayedMatrixStats_1.10.0
[45] ggbeeswarm_0.6.0          viridis_0.5.1
[47] R6_2.4.1                  igraph_1.2.5
[49] compiler_4.0.0


Problem with SPAdes installation - no bin folder I have been trying to download SPAdes for an important project I wanted to do, however, I have encountered a small problem while doing so. It appears that there is no 'bin' folder in my SPAdes installation file.

I have used 7-Zip as my main method of unpacking the file (As I am using Windows, I had to download the actual source code rather than a nice executable), and it came in the .tar.gz format, which I believe 7-Zip is rather good with.

A few pictures detailing my problems are shown below.

I hope I can get a little help on this matter. It is slightly time-sensitive, so I would really appreciate a little bit of assistance shortly. Thank you in advance!

As seen above, the first image is the command I used, the usual test command, the second is the file structure I got (Notice the absence of the 'bin' file), and the third is the error I received from running the spades.py --test command.

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


Retrieval of genomic position by using biomaRt package I have a list of several protein names, their primary gene name and their amino acid sequence (as extracted from the 'SEQRES' section in their corresponding PDB files).

I'm looking for the genomic position (triad of translated nucleotides) of a specific amino acid in the sequence (i.e. residue number 64 in a protein sequence of 344 amino acids).

Is it possible to BLAT an amino acid sequence by using biomaRt in R?
Do you suggest other methods to achieve this task?

BiomaRt cannot do that I'm afraid. There is an Ensembl REST API endpoint that will get you the genomic location of protein coordinates. It needs an Ensembl peptide ID an input though, so you could use the xref endpoint first to get that. There's a bit of sample code in R on those pages, you could use that to script together if you've got a list of these.


Stacked barplot for single cell analysis I want to create a stacked bar-graph with different cell cycles for each cell type within each condition. I have uploaded the file for it. Don't know how to go about it as now I have another condition (cell cycle). I was able to do the bar graph comparing cell cycle per condition but now I want to compare cell type and cell cycle in each condition

output of dput(head(all.combined@metadata))

structure(list(orig.ident = c("treated", "treated", "treated",
    "treated", "treated", "treated"), nCount_RNA = c(1892, 307, 1348,
    3699, 4205, 4468), nFeature_RNA = c(960L, 243L, 765L, 1612L,
    1341L, 1644L), percent.mt = c(0.211416490486258, 1.62866449511401,
    4.45103857566766, 4.4065963773993, 0.0713436385255648, 3.87197851387645
    ), RNA_snn_res.0.5 = structure(c(11L, 11L, 5L, 6L, 11L, 13L), .Label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
    "13", "14", "15", "16", "17", "18", "19"), class = "factor"), seurat_clusters = structure(c(11L, 11L, 5L, 6L, 11L, 13L), .Label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19"), class = "factor"), S.Score = c(0.476893835992198, -0.0200784617568548, -0.0335915198305002, -0.0247184276246385, 0.010785196602457, 0.0190008903712199), G2M.Score = c(0.204441469200986, 0.173804859670862, -0.0313235510969097, -0.0376796363661889, -0.0559526905696905, -0.0122031631356698), Phase = structure(c(3L, 2L, 1L, 1L, 3L, 3L), .Label = c("G1", "G2M", "S"), class = "factor"), old.ident = structure(c(7L,7L, 1L, 4L, 7L, 9L), .Label = c("Fibroblast", "T cell", "Macrophage", "Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell", "Ductal", "Islets of Langerhans"), class = "factor")), row.names = c("treated_AAACGCTAGCGGGTTA-1", "treated_AAAGGTAAGTACAGAT-1", "treated_AAAGTGAGTTTGATCG-1", "treated_AAATGGACAAAGTGTA-1",
    "treated_AACAAAGGTCGACTTA-1", "treated_AACAGGGTCCTAGCCT-1"), class = "data.frame")

output of dput(tail(all.combined@metadata))

structure(list(orig.ident = c("untreated", "untreated", "untreated",
"untreated", "untreated", "untreated"), nCount_RNA = c(901, 823,
1184, 1835, 1147, 1407), nFeature_RNA = c(482L, 479L, 649L, 1043L,
604L, 709L), percent.mt = c(1.77580466148724, 2.91616038882138,
4.22297297297297, 3.86920980926431, 2.0052310374891, 4.05117270788913
), RNA_snn_res.0.5 = structure(c(7L, 7L, 7L, 14L, 7L, 7L), .Label = c("0",
"1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
"13", "14", "15", "16", "17", "18", "19"), class = "factor"),
    seurat_clusters = structure(c(7L, 7L, 7L, 14L, 7L, 7L), .Label = c("0",
    "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",
    "12", "13", "14", "15", "16", "17", "18", "19"), class = "factor"),
    S.Score = c(-0.0320858200243315, 0.0304725660342869, 0.0215996091745327,
    0.0384166213301423, 0.144956251122548, -0.0242770509986111
    ), G2M.Score = c(0.0904224391544142, 0.050148242050667, -0.0178041670730754,
    -0.0112596867977946, -0.0519554524339088, -0.0136533184257381
    ), Phase = structure(c(2L, 2L, 3L, 3L, 3L, 1L), .Label = c("G1",
    "G2M", "S"), class = "factor"), old.ident = structure(c(5L,
    5L, 5L, 5L, 5L, 5L), .Label = c("Fibroblast", "T cell", "Macrophage",
    "Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell",
    "Ductal", "Islets of Langerhans"), class = "factor")), row.names = c("untreated_TTTGGTTGTCTAATCG-18",
"untreated_TTTGGTTTCCCGAGGT-18", "untreated_TTTGTTGAGAACTGAT-18",
"untreated_TTTGTTGAGCTCGGCT-18", "untreated_TTTGTTGAGTGCCTCG-18",
"untreated_TTTGTTGCACGGTGCT-18"), class = "data.frame")

Been using this code to generate the previous graph.

ggplot(CC, aes(x = Condition, y = Percent, fill = Cell_Cycle))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))

Based on your example, you can use count for each phase of cell cycle for each cell types in each condition. Here, I'm using dplyr package to do that but most likely, you could have the same output using various methods:
DF <- rbind(dfhead, dftail)

library(dplyr)
DF_Count <- DF %>%group_by(orig.ident,Phase,old.ident) %>%
  count() %>%
  ungroup() %>%
  group_by(orig.ident,old.ident) %>%
  mutate(Freq = n/sum(n)*100)

# A tibble: 8 x 5
# Groups:   orig.ident, old.ident [5]
  orig.ident Phase old.ident      n  Freq
  <chr>      <fct> <fct>      <int> <dbl>
1 treated    G1    Fibroblast     1 100
2 treated    G1    Stellate       1 100
3 treated    G2M   Tumor          1  33.3
4 treated    S     Tumor          2  66.7
5 treated    S     Mast cell      1 100
6 untreated  G1    Acinar         1  16.7
7 untreated  G2M   Acinar         2  33.3
8 untreated  S     Acinar         3  50

As you can see, DF_Count has the frequency for each cell types of each phase of the cell cycle in function of the condition. We can use DF_Count to get the following plot. Using facet_wrap, you can create two panels based on the condition column and thus represent the cell cycle of each cell types in function of the treatment condition:
ggplot(DF_Count, aes(x = old.ident, y = Freq, fill = Phase))+
  geom_col()+
  geom_text(aes(label = paste(round(Freq, 2),"%")),position = position_stack(vjust = 0.5))+
  facet_wrap(~orig.ident)

Does it look what you are expecting ?
NB: DF here is the addition of your head and tail of your all.combined dataframe.


multiBamSummary output chromosome region annotate issue I was looking how to take readcount from ATAC seq data , this post which suggested to use bedtools and featurecounts to take out readcounts .But multiBamSummary seems to be pretty straight forward. Now I have my count file but only issue is there are no gene names. So now is it possible to annotate these chromosome regions? if so how ?

Any suggestion or help would be highly appreciated

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


IMPUTE2 missing filter option I went through the manuals and can't find the meaning of the below flag:

-filt_rules_l 'filter==0'

Here is the link to manual for this flag: filt_rules_l

Can anyone explain?

IMPUTE version:

impute_v2.3.1_x86_64_dynamic

Edit: The full IMPUTE2 code is irrelevant, I am not using the command to impute, I received an imputed data and the script used to impute. This one flag was a mystery to me. I will update the post with impute software version when I find out.

I have never used this program, but reading the manual in the link it has to mean the following.

The manual says:

  Writing -filt_rules_l 'eur.maf<0.05' on the command line would tell the program to remove any variants with eur.maf values less than 0.05 from the reference panel.

So, -filt_rules_l 'filter==0' would mean that the program removes any variants with filter values of 0 from the reference panel. In the legend file there should be a column named filter, which are filtered out when the value is 0.

Here is the example legend file, in your case there must be a custom additional column named "filter":

$ head ALL_1000G_phase1integrated_v3_chr1_impute.legend

id position a0 a1 afr.aaf amr.aaf asn.aaf eur.aaf afr.maf amr.maf asn.maf eur.maf
rs58108140 10583 G A 0.0407 0.1713 0.1311 0.2071 0.0407 0.1713 0.1311 0.2071
rs189107123 10611 C G 0.0142 0.0276 0.0140 0.0211 0.0142 0.0276 0.0140 0.0211


samtools view: writing to standard output failed: Broken pipe I have the following alignment script which uses BWA:

#!/bin/bash
#usage: sh align_pbs.sh /scratch/waterhouse_team/benth/dbg2olc-40x/good/ output.racon2.fasta bam-soap

mkdir $3

for r1 in $(find $1 -name "*R1*.gz");
do
  base=${r1%_R1*}
  output=$(basename $(echo $r1 | sed 's/R1//g'))
  r2=$(echo $r1 | sed 's/R1/R2/g')

  #Get read group infomration:
  #Source: https://www.biostars.org/p/280837/#310132
  header=$(zcat $r1 | head -n 1)
  id=$(echo $header | head -n 1 | cut -f 1-4 -d":" | sed 's/@//' | sed 's/:/_/g')
  sm=$(echo $header | head -n 1 | grep -Eo "[ATGCN]+$")
  #echo "Read Group @RG\tID:$id\tSM:$id"_"$sm\tLB:$id"_"$sm\tPL:ILLUMINA"
  echo "@RG\tID:${output}\tSM:${output}\tLB:${output}\tPL:ILLUMINA"

  #cat <<EOF
  qsub <<EOF
#!/bin/bash -l

#PBS -N $output
#PBS -l walltime=48:00:00
#PBS -j oe
#PBS -l mem=70G
#PBS -l ncpus=12
#PBS -M m.lorenc@qut.edu.au
##PBS -m bea

cd \$PBS_O_WORKDIR

conda activate bwa-hic
bwa mem \
$2 $r1 $r2 \
-t 12 \
-R "$(echo "@RG\tID:${output}\tSM:${output}\tLB:${output}\tPL:ILLUMINA")" | samblaster -r | samtools view -@ 12 -h -f 0x2 -F 2316 - | samtools fixmate - - | sed 's/^/LP1-/' | samtools sort -@ 12 - -o ${3}/${output}.sorted.dedup.bam

EOF

done

Unfortunatelye, this pipeline cuased the following erros:

samblaster: Version 0.1.24
samblaster: Inputting from stdin
samblaster: Outputting to stdout
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1200000 sequences (120000000 bp)...
[M::process] read 1200000 sequences (120000000 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 316829, 6, 4)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (432, 474, 509)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (278, 663)
[M::mem_pestat] mean and std.dev: (471.59, 53.66)
[M::mem_pestat] low and high boundaries for proper pairs: (201, 740)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 1200000 reads in 615.287 CPU sec, 58.490 real sec
samblaster: Loaded 19 header sequence entries.
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
sed: couldn't write 567 items to stdout: Broken pipe
[E::bgzf_flush] File write failed (wrong size)
samtools fixmate: Couldn't write to output file: Broken pipe
[E::bgzf_close] File write failed
[bam_mating] error while closing output file
samtools view: writing to standard output failed: Broken pipe
samtools view: error closing standard output: -1
samblaster: Unable to write to output file.

Thank you in advance,

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

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

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

What were you trying to achieve with this sed command?


Comparing the gene expression data I'm trying to compare the gene expression data reported in the studies documented in GEO for a specific gene expressed in a tissue of my interest in Homo sapiens.
I compared the values reported in 4 studies, in two studies the expression value of gene A lies between 15-30, in the other two studies the range lies between 1-4.

I would like to ask for advice on whether it is meaningful to compare the data from different studies. What are the factors that one should keep in mind while making such comparisons?

Any help would be highly appreciated

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

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


RCurl issue when installing SCnorm I am working on the Rivanna cluster that is of SLURM-type. Because of that I do not have write permissions and need to install the newest R version myself locally, and use it. For it to be successfully built from source, I needed to also build from source libcurl4-openssl-dev, so I downloaded tar.gz file from here:

https://packages.debian.org/ru/jessie/libcurl4-openssl-dev

and did:

mkdir build
./configure --prefix=/scratch/nv4e/curl/build
make
make install

Then, in the same way I built from source R 3.5.1, the latest current version. It was compiled without any errors, but with some warnings. Then, I submitted the R script that uses SCnorm package:

https://bioconductor.org/packages/release/bioc/html/SCnorm.html

Inside of the script I am having:

source("https://bioconductor.org/biocLite.R")
biocLite("SCnorm")
biocLite("SingleCellExperiment")
library(SCnorm)
library(SingleCellExperiment)

Of course, I am submitting it through .slurm file where I have the main command:

Rscript --vanilla scnorm.R

It fails, apparently because SCnorm requires RCurl, but it can not get one, even though in the log file it is trying to access it correctly:

The libcurl library was built from source in the /scratch/nv4e/curl/build folder. The only suspicion that I have here is that maybe its version is wrong and SCnorm needs another libcurl version? Any suggestions would be greatly appreciated.

  Update

I tried reinstalling libcurl library from the official website, using the latest version:

https://curl.haxx.se/

But it did not fix the issue. I have also tried using libcurl from Anaconda but it gave the same error. So, neither 7.61.0 (latest), nor 7.29.0 (Anaconda), nor 7.38.0 (Debian link) are working.

This is a bug affecting the curl compile, it's been reported here and also here, I suggest trying to compile from the newer release of the libcurl source rather than the .deb you linked.


Comparing aligned amino acids to codon I have a set of 4 amino acids which i have aligned,but i want to compare then with respective  nucleotide sequence in terms of change codon level.
I have the alignment file where i want to take only for the C-terminal domain region region only the CTD starts with amino acids “NITNLC” till it ends with “HAPATV” and rest of the sequence i dont want to take.
In short i want to take the part sequence which starts from "NITNLC" to "HAPATV" then i would like to compare what are the changes in the nucleotide sequence at codon level.

To start with I have used the seqinr library to generate nucleotide alignment from aligned protein sequence.
reverse.align {seqinr}

Then the next step that i want to compare each of my amino acid sequence to my reference sequence which in my case is P0DTC2to rest of the sequence for example I want to compare P0DTC2 to K9N5Q8to the above mentioned CTD region which starts from NITNLC and ends with HAPATVand find what are the codon level changes if there is any I would like to report both in amino acid level and codon level.

First part I guess my approach is correct. Second part I'm not sure how to proceed i guess it's much more than simple parsing I guess!!.
Any help or suggestion would be really appreciated ,and if there exist a R based solution that would be really welcome.
Files
Amino acid file nucleotide file Aligned amino acid file Reverse alignment nucleotide file using seqinr

To set the ball rolling, I tried to implement it manually with your files. Overall, there are three steps:

Locate NITNLC (or HAPATV) in protein P0DTC2.
Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1.
Print the amino acid and the DNA codon of both proteins from step 2.

It works, but only for the first 60 amino acids. I wonder if it's due to the amino sequences in AMINOOO_seq_removed.fasta repeat itself every 60 acids. But why?
#For example, the first three lines of protein P0DTC2
>P0DTC2
MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS
NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV
NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE
...

Step 0: Read the files.
library('seqinr')

#align file containing protein sequences
count_added <- read.alignment('count_added_.clustal_num', format='clustal')
names(count_added$seq) <- count_added$nam

#DNA sequences
rev3.aln <- read.alignment('rev3.aln', format='fasta')
names(rev3.aln$seq) <- rev3.aln$nam

Step 1: Locate NITNLC (or HAPATV) in protein P0DTC2. Again, there are two repeats of NITNLC 60 acids apart (at 807 and 870).
locate <- function(seq, find)
  {address <- gregexpr(paste(strsplit(find, '')[[1]], collapse='[^a-z]*'), seq)
  #substr(seq, address[[1]][1], address[[1]][1]+attr(address[[1]], 'match.length')[1]-1)

  return(list(start=as.numeric(address[[1]]),
             end=as.numeric(address[[1]] + attr(address[[1]], 'match.length') - 1)))
  }

locate(seq=count_added$seq[['P0DTC2']], find='nitnlc')
#start
#807 870
#end
#812 875
locate(seq=count_added$seq[['P0DTC2']], find='hapatv')
#start
#1212 1279
#end
#1217 1285

Step 2: Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1. Instead of using the range 812 - 1279, I chose range 1 - 20 for demonstration purposes.
compare <- function(seq1, seq2, after=0, before=100000)
  {seq1_ = strsplit(seq1, '')[[1]]
  seq2_ = strsplit(seq2, '')[[1]]
  ind = which(seq1_ != seq2_ & grepl('[a-z]',seq1_) & grepl('[a-z]',seq2_))
  ind = ind[ind>after & ind<before]
  #seq1_[ind[1]]
  #seq2_[ind[1]]

  return(ind)
  }
compare(seq1=count_added$seq[['P0DTC2']], seq2=count_added$seq[['K9N5Q8']], after=1, before=20)
# [1] 5  7  8  9 10 13 14

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34

Step 3: Print the amino acid and the DNA codon of both proteins from step 2. Note that ind is based on the align file.
print_amino_codon <- function(ind, seq, seq_gene)
  {locate_amino <- gregexpr('[a-z]', seq)[[1]]
  if (!ind %in% locate_amino) return(NA)
  ind2 = match(ind, locate_amino)

  return(c(amino=substr(seq, ind, ind), codon=substr(seq_gene, ind2*3-2, ind2*3)))
  }

codon(ind=5, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "v" "atg"
codon(ind=6, vseq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "f" "ttc"
codon(ind=7, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "l" "ttg"

codon(ind=5, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "m" "atg"
codon(ind=6, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "f" "ttt"
codon(ind=7, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "v" "gtt"

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34
#<print>         ^^^

#K9N5Q8 gene codon
#gtg ttt cta ctg atg ttc ttg tta aca
#                ^^^ ^^^ ^^^
#P0DTC2 gene codon
#atg ttt gtt ttt ctt
#^^^ ^^^ ^^^


Converting FASTA Sequences to SNP FASTA Format I have two FASTA sequences. Each is from a different species' mitochondrial genome. The sequences are pasted below.
>speciesA_control_region_partial_sequence_mitochondrial
TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA
>speciesB_control_region_partial_sequence_mitochondrial
TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA

My goal is to extract all SNPs in these sequences and format them in SNP FASTA format, as shown below:
>rs11540652
GAACC[A/G]GAGGCCC

I specifically need this format (rather than, say, VCF) so I can search for RFLP sites using the following website:
http://bio.kuas.edu.tw/snp-rflping2/rflpUI.jsp
Does anyone know of an existing program which will accomplish this? Also, what additional materials do I need to detect these SNPs? I presume that I need assembled mitochondrial reference genomes for each species - anything else?

Your question is not very clear. Where does the flanking sequence in :
>rs11540652
GAACC[A/G]GAGGCCC

come from? I'll assume its meant to be the reference sequence and that speciesA is the reference and that the first nuc in the [] is from the reference... I don't know an existing existing program but here's a bespoke one.
speciesA='TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
speciesB='TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
counter=0
for i, nuc in enumerate(speciesA):
    if nuc != speciesB[i]:
        counter+=1
        print ('>rs'+str(counter))
        print (f'{speciesA[i-5:i]}[{nuc}/{speciesB[i]}]{speciesA[i+1:i+6]}')

Gives:
>rs1
AAAGT[G/A]ATGTA
>rs2
TGTAA[T/C]CCTAC
>rs3
CAAAA[G/-]TTATA
>rs4
TTATA[A/G]GAACT
>rs5
CATTA[T/C]ATCAT
>rs6
TTATA[T/C]CATCA
>rs7
AACCA[C/T]TAAAC
>rs8
CCACT[A/G]AACTT
>rs9
CTAAA[C/T]TTCAA
>rs10
AACTT[C/T]AAACT
>rs11
TCAAA[C/T]TCCAT
...

As for 'I presume that I need assembled mitochondrial reference genomes for each species - anything else?' These sequences have clearly been aligned. So there's that.
I hope that helps.


How can I use UpSet plot for expression/RNA-seq data I have a data.frame of 27 000 x 17 and I want to feed it into UpSetR program, but getting the error:

  Error in start_col:end_col : argument of length 0

My data.frame looks like that:

And the upSetR code that I am using:

upset(df, nsets = 4, number.angles = 30, point.size = 3.5, line.size = 2,
mainbar.y.label = "Gene Intersections", sets.x.label = "Genes Per Group",
text.scale = c(1.3, 1.3, 1, 1, 2, 0.75))

What am I doing wrong here? Any suggestions would be greatly appreciated.

  Update

I found here a similar issue:

https://plus.google.com/109635504520741246387/posts/Etujrp6sU1Q

They say that the input should be binary. But then, how can I calculate this binary matrix from the rlog-data? I am not sure...

You need to generate a list of gene ID vectors to feed into upset. In the example below, I am selecting observations (genes) which have an level above the median level. You could also choose only genes above the 3rd quartile, etc. You then need to convert these from a list of character vectors to a presence/absence matrix. UpSetR has a fromList function that does just this.

In future, be more specific about what you want to achieve and supply example data for people to work with.

library("UpSetR")

tmp <- mtcars

tmp$car <- rownames(mtcars)
list <- lapply(colnames(mtcars), function(x) {
  tmp[tmp[[x]] > median(tmp[[x]]), "car"]
})
names(list) <- colnames(mtcars)

upset(fromList(list))

Created on 2019-10-17 by the reprex package (v0.3.0)


Searching for protein motifs with MEME I am trying to look for motifs in a set of proteins returned by a mass spec experiment using MEME. The experiment was carried out in the yeast Candida albicans.

Does anyone know if it is worth calculating either a custom background model or a position specific prior using the non-hit proteins? I don't generally bother for nucleotide motif searching, but I've never done protein motif searching before.

Also advice of whether to use the ZOOPS or ANR models appreciated: is there any reason not to use ANR other than run time?

Some answers by the author Charles Grant via the MEME Suite Google group may be useful to you.

Via: https://groups.google.com/forum/#!topic/meme-suite/yNascbE8Tig/discussion

  Providing your own higher-order background model can greatly improve MEME's ability to discover motifs...

  The command-line utility fasta-get-markov, included in the MEME Suite download, is used to generate custom background Markov models. The input to fasta-get-markov is a FASTA file containing "background" sequences. Ideally, these "background" sequences will be different from the sequences you are analyzing with MEME, but as similar in nature as possible. For example, if you wanted to discover motifs in certain intergenic regions, you might use other sequence data from other intergenic regions to generate the background. The larger the set of "background" sequences is, the better the results will be...

  Typically, you should not specify an order larger than 3 for DNA sequences, or larger than 2 for protein sequences. However, if your input sequences contain higher-order non-random effects that are getting in the way of motif finding, you can follow the following "rules of thumb":

  Use a background model at least four orders less than the shortest motifs you are looking for. So, if you want to find motifs as short as six, I wouldn't use a model higher than order two.
  For an accurate model of order N, you need to use a FASTA file as input to fasta-get-markov with at least 10 times 4**(N+1) DNA characters in it.

  So,

  order-3 requires 2560 characters
  order 4 requires 10240 characters
  order 5 requires 40960 characters etc.

I'm not certain if that input requirement scales differently for proteins (i.e., 20^(N+1)). That may be worth following up on the Google group.

And via: https://groups.google.com/forum/#!topic/meme-suite/UiDhdY6WlK4/discussion

  You should choose the motif occurrence model based on your understanding of the experimental protocol that generated the sequence data and the underlying biology. MEME is trying to find short sequences that are statistically over-represented in your sequence data. To do this, it has to assume a model for how many occurrences of a motif there will be in each sequence. The nature of your experiment should be the basis for the model you choose. For example, data from a ChIP-Seq experiment will typically contain exactly one occurrence of a motif per sequence. On the other hand, sequences drawn from a database of upstream regions of co-regulated genes might contain any number of motif occurrences, including 0.

Via: https://groups.google.com/d/topic/meme-suite/QJiJsy1QxYk/discussion

  The MEME algorithm run time is cubic with respect to the number of input sequences, therefore, it is unsuitable for OOPS (only one per sequence) analyses that have more than 1000 sequences. For ZOOPS (zero or one per sequence) and ANR (any number of repeats) analyses, having more than 1000 sites are also intractable in terms of run time. For example, on a 3 GHz CPU, initialization of the p-value table takes about 10 min for 1000 sites/sequences in OOPS, about a day for 5000 sites, and about 1 week for 10,000 sites. Then the quadratic run time with respect to the total number of characters kicks in next.

Thus it seems like the ZOOPS and ANR models have similar runtime complexity, and use of the MPI-based parallel meme_p binary may be helpful, depending on your analysis:

  The parallel version of MEME scales up to about 128 processors. Please see http://www.sdsc.edu/~tbailey/papers/cabios96.pdf for a discussion of the parallel program. You must bear in mind however, that doubling the total number of sequences input to MEME means that you will need 8 times more processors for the job to finish in the same amount of time.


Which correlation method to compute the correlation score between different clusters of Sc-RNAseq data? I'm analyzing single cell rna-seq data and trying to compute the correlation score between different clusters. Wondering how to choose the correlation method("pearson" (default), "kendall", or "spearman") of the Cor function?

I'm a new comer in bioinformatics and know little about these arithmetics. Can anyone give me some advice or show me the code to decide which method to choose for my data matrix?

Below is my gene expression matrix: rows are gene names,and columns are clusters.I use the default Peasrson method to do analysis,but I really don't know if it is right.

CorOb.cor.exp <- as.data.frame(cor(CorOb.av.exp))
CorOb.cor.exp$x <- rownames(CorOb.cor.exp)
neupeptide.cor.df <- tidyr::gather(data = CorOb.cor.exp, y, correlation,CorOb.cor.exp$x)
a <- ggplot(neupeptide.cor.df, aes(x, y, fill = correlation)) + geom_raster()+scale_fill_gradientn(colors = c("blue", "white", "red"))+
  theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                   size = 12, hjust = 1))+
  coord_fixed()

a

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

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

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

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


Multiway intersection of 8 lists of genes I have 8 lists of differentially expressed genes for 8 time points (2, 4, 6, 8, 10, 12, 14, 16 hours), the number of genes in each of these lists are not equal.

For example there are 630 common genes between h2 and h4 time points but I have 1143 genes in h2 list and 768 genes in h4 list , so common genes in h2 vs h4 would be 630/1143 and common genes h4 vs h2 would be 630/768.

That would would be great if I have a code to produce a matrix or plot of all possible pairwise combinations of my 8 time points.

EDIT: Sorry, I did not see the R tag and wrote the code in bash. I still left the bash example below and added an R version.

Depends on what you want to do, but there are two common ways in which you could store this. Either as wide format:

      h2  h4  h6 ...
Gene1 1   1   0
Gene2 0   1   1
Gene3 1   1   1
...

Where 1 denotes presence and 0 absence, or as long format:

      timepoint
Gene1 h2
Gene1 h4
Gene2 h4
Gene2 h6
Gene3 h2
Gene3 h4
Gene3 h6
...

The best format probably depends on what you want to do. Since you did not specify a method/language, I took the liberty to write a few lines of bash to go from the gene lists to the matrices.
Bash examples assume your files are named h2, h4 etc and are tab separated with gene names in the first column. R examples assume tables is a list of your data frames and timepoints is a vector c("h2", "h4", ...).

Generate table in wide format:

Concatenate gene names of all tables to a separate file
Extract unique values in this file
Perform outer join between the separate file and each table

R:

all_names <- c()
for (tp in timepoints){
  all_names <- append(all_names, tables[[tp]]$genes)
}
all_names <- unique(all_names)

wide_table <- matrix(nrow=length(all_names), ncol=(1+length(timepoints)))
wide_table <- as.data.frame(wide_table)
colnames(wide_table) <- append("genes", timepoints)

for (tp in timepoints){
 tmp <- tables[[tp]]
 tmp$tp <- 1
 tmp <- merge(wide_table, tmp, by="genes", all=T)
 wide_table$tp <- tmp$tp
}

Bash:

Note: you will might need to add headers to this table

for timefile in $yourfiles;do
  # make a file with all gene names
  awk -v tp=$timefile '{print $1}' >> names_file
done

sort names_file | uniq > tmp && mv tmp names_file

for timefile in $yourfiles;do
  join -a1 -a2  -j 1 -o 0,1.2,2.2 -e 0 $names_file <(awk '{print$1,1}' $timefile | sort)
done

Generate table in long format:

For each table, add a column with timepoint
Concatenate all tables

R:

df <- data.frame()
for (tp in timepoints){
  tables[[tp]][[timepoint]] <- tp
  df <- rbind(df, tables[[tp]] )
}

Bash:

   for timefile in $yourfiles;do
      awk -v tp=$timefile '{print $1,timefile}' >> long_table
    done

EDIT2: To follow up on this, here is a simple way (using the package dplyr) to count the number of common genes between all possible combination of time points. Assuming you have your data in long format:

df <- data.frame(genes=c("G1","G1","G2","G3","G3"), tp=c("h2","h4","h2","h4","h6"))

  genes tp
1    G1 h2
2    G1 h4
3    G2 h2
4    G3 h4
5    G3 h6

You can get the co-occurrences using :

library(dplyr)
# Make timepoint ordinal factor
df <- mutate(df, tp = factor(tp, ordered=T))
df <- df %>%
  # Compute all pairwise co-occurrences
  full_join(df, by="genes") %>%
  group_by(tp.x, tp.y) %>%
  summarise(N_genes = length(unique(genes))) %>%
  # Exclude self pairs (gene in h2 and h2)
  filter(tp.x != tp.y) %>%
  # Check for duplicate pairs in reverse order (h2-h4 and h4-h2)
  mutate(comb=ifelse(tp.x < tp.y,
                     yes = paste(tp.x, tp.y),
                     no  = paste(tp.y, tp.x)))
# Remove the duplicate pairs
df <- df[which(!duplicated(df$comb)),]

And visualise it using for example a heatmap:
(ofc it looks stupid because I have a tiny sample where not all combinations are even present)

# Visualize number of genes present in each combination
library(ggplot2)
ggplot(df, aes(x=tp.x, y=tp.y, fill=N_genes)) + geom_tile()


Selecting sites from VCF which have an alt AD > 10 I have high-depth variant calling created using the HaplotypeCaller with --output_mode EMIT_ALL_SITES I'm interested in finding all sites (regardless of genotype call heterozygous or homozygous) where at least one of the alternative alleles have an AD value (Allelic Depth) greater than 10, I.e. are supported by more than 10 reads. Also ideally I want back more than just the first alternative allele.  Note that I don't want back lines of VCF were we only see an AD count for the ref allele only.
So in the example VCF snippet below I'm wanting to select lines: 6,7,8,12,13 and 14, which have GT:AD values 1/1:1,988:989 0/1:116,92 0/1:220,234 0/1:62,611 1/1:0,109 respectively.
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  12908_DIAG
3   187446740   .   T   .   Infinity    .   AN=2;DP=1095;MQ=60.00   GT:AD:DP    0/0:1095:1095
3   187446741   .   C   .   Infinity    .   AN=2;DP=1117;MQ=60.00   GT:AD:DP    0/0:1117:1117
3   187446752   .   A   .   Infinity    .   AN=2;DP=1297;MQ=60.00   GT:AD:DP    0/0:1297:1297
3   187446763   .   C   .   Infinity    .   AN=2;DP=1494;MQ=60.00   GT:AD:DP    0/0:1494:1494
3   187451574   .   C   .   Infinity    .   AN=2;DP=1493;MQ=60.00   GT:AD:DP    0/0:1493:1493
3   187451609   rs1880101   A   G   39794.03    .   AC=2;AF=1.00;AN=2;BaseQRankSum=1.859;ClippingRankSum=0.000;DB;DP=995;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=0.000;QD=24.56;ReadPosRankSum=0.406;SOR=8.234  GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0
4   1803279 .   T   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0.000;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0.00;MQ=59.97;MQRankSum=0.000;ReadPosRankSum=-2.523;SOR=6.357   GT:AD:DP:GQ:PL  0/0:211,23:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0.000;DB;DP=215;ExcessHet=3.0103;FS=1.110;MLEAC=1;MLEAF=0.500;MQ=59.97;MQRankSum=0.000;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809  GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803671 .   C   A   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-0.880;ClippingRankSum=0.000;DP=450;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.953;SOR=0.572    GT:AD:DP:GQ:PL  0/0:445,2:447:99:0,1272,15958
4   1803681 .   T   C   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.654;ClippingRankSum=0.000;DP=483;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.422;SOR=0.664    GT:AD:DP:GQ:PL  0/0:479,2:481:99:0,1408,18538
4   1803703 .   A   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.704;ClippingRankSum=0.000;DP=458;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=0.299;SOR=0.497 GT:AD:DP:GQ:PL  0/0:454,2:456:99:0,1325,18095
4   1803704 rs2234909   T   C   6676.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0.000;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=14.71;ReadPosRankSum=0.324;SOR=0.849   GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=8.083;ClippingRankSum=0.000;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=16.51;ReadPosRankSum=0.180;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
4   1805296 rs3135883   G   A   3876.03 .   AC=2;AF=1.00;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.22;SOR=9.401    GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0

A dropbox link for file
I'd initially considered using GATK's SelectVariants but I'm not sure JEXL has the ability to select out what I want specifically other than a blanket AD > 10 which will give me both ref and alt alleles with AD > 10.  Perhaps there is a bioawk solution or something more elaborate with coreutils which could successfully return sites with an alt AD count > 10?

using vcfilterjs

and the following script:

function accept(vc)
    {
    var i,j;
    for(i=0;i< vc.getNSamples();++i)
        {
        var genotype = vc.getGenotype(i);
        if(!genotype.hasAD()) continue;
        var ad = genotype.getAD();
        /* loop over AD starting from '1' =  first ALT */
        for(j=1  ;j< ad.length ;++j)
            {
            if(ad[j]>10) return true;
            }
        }
    return false;
    }
accept(variant);

usage:

java -jar dist/vcffilterjs.jar -f script.js Test.vcf


How do I compare measured ratios to calculated ratios in peptide mass spectrometry? I have two different sets of internal standards in a sample (1), and a sample (2) with these internal standards as well as about 30 peptides. Each peptide has one of each of these kind of standards. From sample 1, I can get the true ratio between standard 1 (S1) and standard 2 (S2), and from sample 2 I can get the ratio between the analyte (L) to each of the standards:

Sample 1: S1/S2

Sample 2: L/S1, L/S2

I can thus calculate the ratio S1/S2 from the measured ratios in sample 2. Each sample is run in 5 technical replicates (that is; 5 injections into the LC-mass spectrometer), so S1/S2 ratios for one particular analyte could look like this (mock data):

measured    calculated
0.967       0.987
1.007       0.967
1.044       1.012
1.041       1.025
1.048       1.046

I want to compare the measured ratios to the calculated, to see if there are any bias in my acquisition of data. As the calculated and the measured ratio are from different samples, the replicates aren't "linked" so I could not use for example Bland Altman analysis. Should I use some kind of ANOVA? Or do you have any other recommendation?

I would also like to make an overall assessment of the ratios of all 30 peptides to see if the bias is different for different peptides. How could I do that?

Any input would be greatly appreciated.

From a statistical perspective ratios are icky because the variance of the sampling distribution of a ratio decreases with the size of the numerator and denominator being sampled to create the ratio. If you can get at the numerator and denominator separately, then you should absolutely do that and look at count based tests like the Chi-square test, Fisher's exact test, or Poisson regression.

That said, a t-test would be appropriate in the two-sample case if its assumptions are met.  You need to check that both samples are normally distributed, ideally via a qq-plot. You also need to check the variance of both samples and choose a variant of the t-test that does or does not assume equal variances.

If the samples are not normally distributed, then the next approach would be a Mann-Whitney U test.


Correct for gene length or read counts in GO enrichment analysis It is a well reported fact that GO analysis of RNAseq results is affected by a number of biases, including length bias and expression level bias.

The bioconductor package goseq allows you to correct for these biases.

By default it corrects for length bias, but you can also get it to do read count bias. Using read counts to do the correction is attractive because in theory it should account for both sources of bias ($read counts\approx expression \times length$).

I'm doing an enrichment analysis were I have tried both options (length and read counts) and get very different answers. If I run a binomial regression on expression and length vs probability of being differential, I can see that both are independently important.

> model <- glm(sig ~  expression + log(length), data=retained_genes,  family=binomial(link="logit"))
> print(anova(model, test="Chisq"))

Analysis of Deviance Table

Model: binomial, link: logit

Response: sig

Terms added sequentially (first to last)

                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
NULL                                    6676     4507.1
expression              1  114.998      6675     4392.1 < 2.2e-16 ***
log(length)             1  102.553      6674     4289.5 < 2.2e-16 ***
expression:log(length)  1   34.094      6673     4255.4 5.252e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

My understanding of the subject is that the bias of the gene length (and other bias) should be taken care when analyzing the expression and before enrichment analysis. The enrichment analysis should be done once the corrections has been performed. Because as the abstract of the GOseq paper states:

  GO analysis is widely used ... but standard methods give biased results on RNA-seq data due to over-detection of differential expression for long and highly expressed transcripts.

So first, take care of the differential expression bias by length and then use the GO to reduce complexity. How you take care of the bias in the RNA-seq data is another question. But the cqn package of Bioconductor can correct the expression by gene length and GC content. However, this correction might hurt the differential tool used (see this discussion in Bioconductor), so it might be better in some cases to use GOSeq.

Now, to the question itself:

  So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Use whatever correction method that yields better differentially expressed genes (DEG. If you find that the correction for length improves the accuracy of the predictions of DEG better than correcting by length and GC, then use that one.

Another option to obtain accurate GO terms, then you could use other testing procedures which don't rely uniquely in the Fisher test, such as those that take into account the structure of the GO graph. TopGO use this approach (note that it is a bit difficult to work with this package),this will reduce the role of the gene length bias (and probably other bias) in the resulting significant GO.


Does picard markduplicate toggle PCR duplicate samflag I have a RNA-seq bam file and there are few reads that are puzzling me.

According to the bam header, this bam file is sorted by coordinate, created using tophat and markduplicate step is not done. But some reads are marked for being duplicate in the samflag. What is worse is when I run picard markduplicate, these reads' pcr duplicate flag is toggled, marking them not a duplicate. Also I manually found the duplicate of this read (identical reads with same start positions and mates start position) so initial marking looks true.

So my questions are:

Any idea why would this happen?

Does Tophat mark reads for being duplicate? ( I don't think so)
And does picard markduplicate toggle the duplicate flag if the reads are already marked as being duplicates?

Here is how the read looks before and after mark duplicate step.

Before:

C0RTF 1187 17 7579880 255 61M10754N40M = 7579927 10902 CTC...
0UNP1 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...

After Markduplicate:

C0RTF 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...
0UNP1 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...

Yes, if MarkDuplicates encounters a pair that's marked as a duplicate that it considers (for whatever reason) to not be a duplicate then it will unset the duplicate mark. You can test this yourself by making a small BAM file either with or without duplicate entries but some marked regardless. In the case of actual duplicates, the duplicate flag will be maintained. In the case where there are no actual duplicates then the flag will be removed. At least that's the case in version 2.14.


bcl2fastq2 v2.20.0 converstion software on Ubuntu 18.04 bcl2fastq is a program that converts Illumina run output from .bcl to fastq. I'm having difficulties installing the software I downloaded from Illumina on my Ubuntu 18.04. I've followed their installation manual, but the software is not building when I enter the command "make".  I receive the error message... "make: *** No targets specified and no makefile found. Stop. Has anyone had similar experience? Is there another software that will install easier on Ubuntu?

Software:
https://support.illumina.com/downloads/bcl2fastq-conversion-software-v2-20.html

Manual:
https://support.illumina.com/content/dam/illumina-support/documents/documentation/software_documentation/bcl2fastq/bcl2fastq2-v2-20-software-guide-15051736-03.pdf

Steps followed:

export TMP=/tmp
export SOURCE=${TMP}/bcl2fastq
export BUILD=${TMP}/bcl2fastq2-v2.20.0-build
export INSTALL_DIR=/usr/local/bcl2fastq2-v2.20.0

cd ${TMP}
tar -xvzf bcl2fastq2-v2.20.0.499-Source.tar.gz

mkdir ${BUILD}
cd ${BUILD}
chmod ugo+x ${SOURCE}/src/configure
chmod ugo+x ${SOURCE}/src/cmake/bootstrap/installCmake.sh
${SOURCE}/src/configure --prefix=${INSTALL_DIR}

cd ${BUILD}
make

The final command is supposed to be make install, but the configure step fails with:

[ . . . ]
-- Configuring incomplete, errors occurred!
Couldn't configure the project:
/tmp/bcl2fastq2-v2.20.0-build/bootstrap_cmake/bin/cmake -H'/tmp/bcl2fastq/src' -B'/tmp/bcl2fastq2-v2.20.0-build' -G'Unix Makefiles' -DBCL2FASTQ_NAME_SHORT:STRING="bcl2fastq" -DBCL2FASTQ_NAME_LONG:STRING="BCL to FASTQ file converter" -DBCL2FASTQ_COPYRIGHT:STRING="Copyright (c) 2007-2017 Illumina, Inc." -DBCL2FASTQ_VERSION_MAJOR:STRING="2" -DBCL2FASTQ_VERSION_MINOR:STRING="20" -DBCL2FASTQ_VERSION_PATCH:STRING="0" -DBCL2FASTQ_VERSION_BUILD:STRING="422" -DBCL2FASTQ_VERSION:STRING="2.20.0.422" -DBCL2FASTQ_SOURCE_DIR:STRING="/tmp/bcl2fastq/src" -DBCL2FASTQ_PREFIX:PATH="/usr/local/bcl2fastq2-v2.20.0" -DBCL2FASTQ_EXEC_PREFIX:PATH="" -DCMAKE_INSTALL_PREFIX:PATH="/usr/local/bcl2fastq2-v2.20.0" -DBCL2FASTQ_BINDIR:PATH="" -DBCL2FASTQ_LIBDIR:PATH="" -DBCL2FASTQ_LIBEXECDIR:PATH="" -DBCL2FASTQ_INCLUDEDIR:PATH="" -DBCL2FASTQ_DATADIR:PATH="" -DBCL2FASTQ_DOCDIR:PATH="" -DBCL2FASTQ_MANDIR:PATH="" -DCMAKE_BUILD_TYPE:STRING="RelWithDebInfo" -DCMAKE_PARALLEL:STRING="1"
Moving CMakeCache.txt to CMakeCache.txt.removed

Install it with conda.
conda install -c dranew bcl2fastq


What are doublets in single cell RNA-seq data? I am reading The Tabula Muris Consortium et al. (pp).

  In some organs, cells with more than 2 million reads were also excluded as a conservative measure to avoid doublets.

How exactly is a “doublet” defined? For example, is doublet a set of cells sequenced as a single cell (and so, the number of transcripts is double)?

Is doublet a set of cells sequenced as a single cell?

Yes. Depending on the method of single cell sequencing it may be more or less likely for groups of cells to be captured and barcoded with the same "unique" barcode. This is more likely in split-pool RNA sequencing (e.g. SPLiT-seq), and less likely in cell-capture RNA sequencing (e.g. Fluidigm C1). Doublets can also be created through physical / experimental processes (e.g. from tissues that were not completely dissociated).
Bear in mind that detecting doublets is not as simple as counting for doubling of transcript expression, because the expression profiles are different in different cells. That's why single-cell sequencing is useful in the first place.


Measure the purifying selection for certain taxa along a phylogeny I am posting this message because I need clarity about the analysis I want to perform.
In my analysis, I have a homologous gene in 13 species, and I would like to evaluate the selection pressure of this gene (its dN/dS) in 3 species of the tree and test if this value is significantly different from the neutral model (dN/dS = 1).
So here is the protocol I am considering, and I would like you to tell me if I am right or not:
To calculate this selection pressure, the ratio of non-synonymous to synonymous substitution rates (dN/dS) must be calculated.
For this I use CODEML implemented in ete3.
So I need :
A phylogenetic tree of the gene (performed with the AA alignment of the gene via Iqtree)

A codon alignment (performed with the DNA alignment of the gene with pal2nal)

Mark branches in the tree.

Here is the exemple where I want to get the dN/dS value of SP5,6 and 7.:
         /-SP1
        /-|
       |   \-SP2
       |
     /-|      /-SP3
    |  |   /-|
    |  |  |   \-SP4
    |   \-|
    |     |   /-SP5
  /-|      \- |  /-SP6
 |  |         \-|
 |  |           |
 |  |            \-SP7
 |  |
 |  |   /-SP8
-|   \-|
 |      \-SP9
 |
 |     /-SP10
 |  /-|
 | |   \-SP11
 |-|
   |   /-SP12
    \-|
       \-SP13

And I want to calculate the dN/dS of SP5,SP6 and SP7 in order to confirme the hypotesis that the gene is under strong purifying selection on this 3 species.
I mark the SP5, SP6 and SP7, run a free-branch model and look at the Codeml output (the dN/dS value should be the same for the 3 marked species right?) Then this value is kind of the mean dN/dS value between three leafs?
Suppose I get an omega = 0.15
If I want to test that this value is significantly different from 1, I need to mark the same sheets and run the neutral model, then compare the two models (neutral vs. free) using chi2, and if the p-value is <0.5, then the dN/dS value (0.15) is significantly different from 1.
I'm right  ?

The following is not the clearest explanation of dN/dS humankind has seen, but I hope you get the gist and core idea of likelihood is absolutely correct. I am wondering if you are thinking of a different test to the way I would do this calculation. This ain't a text book clarity answer.
Anyway what you are proposing is a likelihood ratio test (LRT) for  selection based on the parameterisation of omega. Your idea isn't quite correct ... its not dN/dS value you are looking at to assess significance, its the likelihood. This point is critical and is the most important difference in our thinking on this test.
LRT works because its a log likelihood and subtracting two log likelihood is the same as dividing them (I forget the theory - but you subtract the likelihood to work out the p-value via Chi squared). There's some theory about why the Chi-squared distribution works but I can't remember it.
You need to work out the degrees of freedom which is the parameter differences between 'free branch model' and the neutral model and this value is LOW. If the neutral model is the "global model", then its 2 degrees of freedom - IN YOUR TREE. This is because the SP you are looking at form a single clade (monophyly). From memory you stick a #1 on the branch leading into the clade, ie. defining the node and thats your 'model' on the treefile. Here you have two dN/dS one for the whole tree excluding your clade and one for your taxa. This gives you a likelihood.
Note you could stick a #1 on each branch within the clade, i.e. the treefile. That means something different - like every branch was subject to ... in this case purifying selection ... not just the common ancestor.
You now work out the likelihood of the 'global' model... i.e. remove the #1 and make sure it gives you only 1 dN/dS value.
Subtract the likelihoods and multiple by two - if its LRT (super important).
You now need to work out the degrees of freedom. If its the same model but one has an additional node with a different dN/dS ... the df = 2. If you've dumped a #1 on every branch in the clade ... I still think df = 2 - you need to check.
If the difference in likelihood results in p < 0.05 you are correct the difference is significant ... whatever that difference might be.
If the neutral model is some different calculation of omega within the population genetics process (see here) then it will be the parameter differences between your model and the 'neutral' model.
To reiterate
I like simulations so much I type in reiterations ... So again the LRT usual will usually subtract one likelihood from another and THEM multiply the probability by two ... and compare the result on the Chi squared table for the difference in degrees of freedom. And again what you are looking at is not the value of dN/dS but the difference in LIKELIHOOD. That is very important. You need to be careful with the degrees of freedom on an LRT because they can introduce an additional subtraction constant, test dependings (its ages since I've done one so I can't remember).
If the test is significant and one dN/dS is x for the clade and y for the rest of the tree, then that differences is significant and if y is >1 (bit unlikely) and x < 1 then sure thats purifying selection. Usually clades evolve via positive selection rather than purifying selection.
From memory that is my understanding ... but if you have a neutral model of omega and a purifying model of omega ... yeah sort of, but usually thats for a 'global' model, i.e. no differences within the tree. In your case your interested in a specific clade. The idea is still the same its measured through likelihood.
Its so long since I've done this, but I do caution not to confuse two different tests.

I do understand Yang et al did an MCMC for dN/dS ... what I'm thinking of here is maximum likelihood. You'd need to watch for that carefully and thats what you are referring to when you say Chi-squared, because that is LRT 100%. LRT is maximum likelihood.


Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern What's the most straightforward way to filter a BAM file for read pairs where one or both of the reads starts with a given sequence pattern? Would a tool that deals with UMIs work for this? Which one would be best?

First option

My solution uses the following steps:

use picard sortsam to sort the records on query-name (not samtools sort because the order is not the same between java and C )
use jjs (java scripting engine) and the htsjdk library to build a bufferof reads having the same name. If any read starts with a pattern, print the whole buffer. The htsjdk should be already contained in the picard.jar.
sort back to coordinate order using samtools.

All in one:

java -jar /path/to/picard.jar SortSam I=in.bam O=/dev/stdout  SO=queryname |\
jjs -cp /path/to/picard.jar  script.js |\
samtools sort -T TMP -o out.bam -

with script.js

var out=java.lang.System.out;
var File = Java.type("java.io.File");
var SamReaderFactory = Java.type("htsjdk.samtools.SamReaderFactory");
var SamInputResource = Java.type("htsjdk.samtools.SamInputResource");
var SAMFileWriterFactory = Java.type("htsjdk.samtools.SAMFileWriterFactory");
var samReader =  SamReaderFactory.makeDefault().open(SamInputResource.of(java.lang.System.in));
var samWriter =  new SAMFileWriterFactory().makeBAMWriter(samReader.getFileHeader(),true,java.lang.System.out);
var iter=samReader.iterator();
var buffer=[];
for(;;)
        {
        var rec=null;
        if(iter.hasNext()) rec=iter.next();
        if(rec==null || (buffer.length>0 && !rec.getReadName().equals(buffer[0].getReadName())))
            {
            var i=0;
            for(i=0;i< buffer.length;++i)
                {
                if(buffer[i].getReadString().match(/^AAAA/)) break;
                }
            if(i!=buffer.length)
                {
                for(i=0;i< buffer.length ;++i)
                    {
                    samWriter.addAlignment(buffer[i]);
                    }
                }
            if(rec==null) break;
            buffer=[];
            }
        buffer.push(rec);
        }

samReader.close();
samWriter.close();

A second option is to use samtools view and awk to export the read names:

 samtools view  input.bam | awk -F '\t' '{if($10 ~ /^AAAA/) print $1;}' | sort | uniq > names.txt

and to use Picard FilterSamReads and READ_LIST_FILE to filter the reads.

UPDATE: 3rd solution:

I've updated my tool samjdk: there is now a new option --pair to work with query-name-sorted BAM files.

The command would be:

java -jar picard.jar SortSam I=S1.bam O=/dev/stdout  SO=queryname |\
java -jar dist/samjdk.jar --samoutputformat BAM --pair \
     -e 'return records.stream().anyMatch(R->R.getReadString().startsWith("AAAAAAG"));' |\
samtools sort -T TMP -o out.bam -


Merging ssu-align alignments ssu-align is an HMM-based MSA tool for ribosome's small subunit (SSU) sequences. It's bundled with three SSU models:

archaeal,
bacterial
eukaryotic

Approach Given a set of archaeal and bacterial sequences the tool will first identify the best model for each sequence and produce two separate alignments for archaea and bacteria, and the tool provides no merging functionality.

Rationale The underlying issue with SSU is the secondary structure, which results in: 1) indels; 2) pairwise point mutations (to maintain the secondary structure)

Problem While you can force it to use a single model for all sequences, aligning archaeal sequences to the bacterial model will produce poor results.

Solution? Nevertheless, the authors of GreenGenes claim that they've used ssu-align to produce two MSAs (one for archaea and bacteria) and build a joint phylogenetic tree using FastTree, but they say nothing about the way they've managed to merge the alignments. Am I missing some functionality from SSU-align?

Note Merging two giant MSAs using conventional tools, such as MAFFT, won't get me far.

P.S.

I've already tried contacting the authors of the GreenGenes paper, but the corresponding author hasn't responded.

I think this is more to do with tree building than alignment construction.

I believe the underlying issue is the 1) distribution of indels, this results in a blockwise deletion in the resulting phylogeny and 2) nucleotide saturation.

Indels Thus, for indels

Archaeal
XXXX----YYYYY
bacterial
----ZZZZYYYYY

The resulting phylogeny will be based on YYYYY and XXXX is ignored.

The alternative explanation is that the genetic distances involved have resulted in saturation of the nucleotide sequences. SSU has to be a nucleotide phylogeny and the level of reversion mutations over a vast amount of evolutionary time can become indistinguishable from random.

With regards the first point, you can replace the indels "-" with "N". This is a cheat to trick a tree building program into not discarding the indel. Model based tree building programs have no ability to create a genetic distance between a base pair and an indel, so this data is essentially discarded from the tree. Worse it will remove the site from phylogenetic inference - which can be a massive loss of data. Two notes one this,

Some max. likelihood, Bayesian tree programs can avoid blockwise deletion, I don't know programmatically how they do this;
In PAML the distribution of indels was modelled, it doesn't appear to have been widely adopted however, but might be of value.

Alternative The alternative method is really for protein genes and called "RY" coding, it might work on SSU because of the compensatory mutations involved in maintaining a secondary structure of the SSU.

Word of warning, the saturation may be so great, very little nucleotide point mutation is preserved. Essentially we are talking about preservation of phylogenetic signal across the origins of life ... and that is a very, very long time ago.


Calculating p-value for introns which are retained (or not) between 2 conditions I am trying to calculate p-values and ultimately FDR values for RI introns for 2 conditions of experiment I have(control and knockdown)

For every condition I have 3 biological replicates each: 3 control replicates and 3 knockdown replicates.

This is the header line of my input file.

chrom   Istart  Iend    strand  avg-L_R_counts-1    Sk_counts-1 avg-L_R_counts-2    Sk_counts-2

-1 and -2 are 2 sample conditions control and knockout respectively

avg-L_R_counts = for a particular intron, this is the average of reads 5bp up stream and downstream of the 5' end and at the 3'end of the intron

Sk_counts = This is the number of reads spanning exon-exon junction for that intron.

I have approximately 350 introns(this is yeast). I can get the Percent Spliced In values but I am interested in getting a p-value and then I can calculate the FDR values.
Can anyone point me as to what approach I should take to calculate a p-value, or what test should I use?

Here is how we do this:

First break your gene model down into non-overlapping chunks, such that any possible gene model can be built from a combination of those chunks. E.g.

transcript 1 |>>>>>>>>>|----------|>>>>>>>>|---|>>>>|
transcript 2 |>>>>>>|-------------|>>>>>>>>|---|>>>>|
transcript 3 |>>>>>>|-------------|>>>>>>>>>>>>>>>>>|
chunks       |+++1++|+2|+++++3++++|+++4++++|+5+|+6++|

this can be achieved using cgat gtf2gtf --method=genes-to-unique-chunks from the cgat-apps toolkit if you don't want to code it yourself.

I call
chunks 1, 4 and 6 constitutive exon
chunk  3 constitutive intron
chunk  2 alternate intron
chunk  5 annotated retained intron.

Now use featureCounts to count the number of reads in each chunk.

Test the chunk counts for usage using DEXSeq.

Filter the results to just include chunks of type 2, 3 and 5 and recompute the FDRs.


Does rRNA depletion protocol give higher number of mapped reads in Intronic regions? Recently, I have downloaded a publicly available dataset, which are 350 tumor samples. I see the following information from the published paper.

They used Ribo Zero Gold and rRNA was depleted. Strand specific data. After aligning the data I did some alignment quality check with Qualimap RNA-Seq QC tool. I visualised the bam files in IGV. Alignment is good. For all samples 90% alignment rate was seen. I observed that in all samples Higher percentage of mapped reads were originating in Intronic regions. Followed by Exonic and intergenic regions.
I have seen a post here Reads mapped to exonic, intronic and intergenic regions where they say high intronic reads could be because of contamination. I googled about higher reads in intronic regions and found some papers Evaluation of two main RNA-seq approaches for gene quantification in clinical RNA sequencing: polyA+ selection versus rRNA depletion and some other links RIBO-DEPLETION IN RNA-SEQ – WHICH RIBOSOMAL RNA DEPLETION METHOD WORKS BEST? in which they said Greater intronic reads were with rRNA depletion protocol.
And even in this RNA-SEQ tutorial, it is mentioned that - A higher intronic mapping rate is expected for rRNA removal compared to polyA selection.
So, my question:
I am working with lncRNAs. So, I'm using the samples prepared with rRNA depletion protocol. Is this higher intronic rate is common in rRNA depleted dataset or do I have to check anything else to proceed further with these samples?

It's not so much that you have "intronic contamination" or "genomic contamination", rather you're not selecting explicitly for full-length mature transcripts with rRNA depletion. That is the most common cause for higher intronic read rates. There's nothing you can do about this post-hoc, just continue along.
BTW, many lncRNA's are polyadenylated, so you'll keep them with poly-A selection.