text:
What is the difference between AnnotationDbi/org.Mm.eg.db and biomaRt/Mus.musculus for converting to gene symbols?
I am interested in the differences between AnnotationDbi/org.Mm.eg.db and biomaRt/Mus.musculus. They yield the same results as seen in the code below.
library(edgeR)
library(biomaRt)
library(Mus.musculus)
library(AnnotationDbi)
library(org.Mm.eg.db)

url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
           "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(i, overwrite=TRUE)

files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
           "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
           "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
           "GSM1545545_JMS9-P8c.txt")

x <- readDGE(files, columns=c(1,3))
geneid <- rownames(x$counts)

genes <- AnnotationDbi::select(org.Mm.eg.db, keys=geneid, columns="SYMBOL", keytype="ENTREZID")
dim(genes)
genes2 <- biomaRt::select(Mus.musculus, keys = geneid, columns = c("SYMBOL"), keytype = "ENTREZID")
dim(genes2)

Do they use different databases? Is one preferred or recommended over the other?

AnnotationDBi uses AnnotationDb objects which is the virtual base class for all annotation packages. It contain a database connection and is meant to be the parent for a set of classes in the Bioconductor annotation packages. Examples of classes derived from the AnnotationDb object include (but are not limited to): ChipDb, OrgDb, GODb, OrthologyDb and ReactomeDb. The database from one of these classes is used to extract annotations. org.Mm.eg.db belongs to OrgDb class.
Mus.musculus  represents a collection of annotation packages that can be used as a single object. It combines different information about same organisms (Mus musculus) into one place. It also uses org.Mm.eg.db as a dependency. So, both uses same database to extract Gene symbols. That is the reason they are producing same results. There is helpful documentation available on Bioconductor.


text:
Effect size in power analysis when dealing with Poisson based variant caller
I was trying to estimate the minimum sequencing depth required to confidently detect somatic variant using Poisson distribution based variant caller (such as somatic variant caller or Pisces from illumina).
Power analysis is used which requires significance level (0.05), power / sensitivity (99%) and effect size to estimate sample size (sequencing depth). I am not sure how to calculate the effect size in such context. Can anyone share some insight? Thanks a lot.

It's standard stats, it's a while since I've done this, but it's textbook stuff

Try running a pilot study to get some preliminary results you can use to estimate the effect size = difference / SD.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Defining methylation level of CpG island
Differentially methylated site is technically measured by each probe. However, CpG island within gene promoter is comprised of numerous probes. Do we have a consensus of how to determine the methylation level of CpG island using methylation level of probes within corresponding CpG island?

In my experience, the most common approach to find DMR is to use the DSS R package to call and merge DML into regions.
DMC/DML (Differential Methylation Cytosine/Loci) is often calculated by comparing the counts of methylated and unmethylated cytosine at one position between samples.
Steps are to first call DMLs with the callDML() function. Then call DMRs from these first results with callDMR(). Finally annotate the DMLs or DMRs using a GFF file. From there you can look for DMR falling into CpG island. Some of the parameters that can be tuned are the minimum region length, minimum number of differential CpG sites, percentage of CpG sites required to be significant. However I have never heard of what would be the standard values for these parameters. It may not be clearly established yet.
So because the process of calling DMLs and DMRs need robust statistical testing, it is suggested to first call the differential methylated states and then annotate it.

Yes, most current studies use DMR to measure methylation level. However, DNA methylation impose its impact by suppressing gene expression. I think it is more reasonable to measure the magnitude of methylation impact on gene expression than merely defining statistically methylated region. To do this, we may need better understanding on the correlation between methylation and mRNA expression (things like whether different CpG site have equal contribution on the expression of corresponding gene)

Well as you may know the effect of methylation on genes expression is far from being clear, most studies showing poor negative correlation between differential gene expression and methylation. An interesting approach was the one of Vanderkraats in https://doi.org/10.1093/nar/gkt482. The main idea was to look for pattern of differential methylation around genes promotor that could explain genes expression.

Thanks @PaulEndymion a detailed adn useful response

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Why did expression based subtypng of breast cancer gain much more acceptance than others
This is may not be entirely technical question but rather a academic question. But the technique behind the application is within the scope of bioinformatics. So I would try to ask here that:
In each cancer type, there have been tons of papers that describe RNA-Seq based subtyping resulting in molecularly distinct subgroups with survival difference. However, only very few of them gained acceptance at clinical setting, such as breast cancer (luminal subtype A, luminal subtype B, normal breast-like subtype, HER-2 overexpression subtype and basal-like subtype).
My question is: why did this molecular subtyping of breast cancer gain much acceptance but not others? Did this one own some advantages that other did not?

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

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

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

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


text:
Error in running CIRIquant package for 2 RNA-seq in Ubuntu 16.04
"CIRIquant" is a Python package for detecting Circular RNA. I ran "CIRIquant" individually and it was run. consequently, I ran CIRIquant through the below command and parameters in the terminal:
 ~/.local/bin/CIRIquant -t  4

     -1 ./SRR1176949_1.fastq.gz
     -2 ./SRR1176949_2.fastq.gz
     --config ./CIRIquant-Config.yml
     -o ./SRR1176949
     -p  SRR1176949

But I got the below message:
Traceback (most recent call last):
  File "/home/genetic/.local/bin/CIRIquant", line 11, in <module>
    sys.exit(main())
  File "/home/genetic/.local/lib/python2.7/site-packages/CIRIquant/main.py", line 89, in main
    config = check_config(check_file(args.config_file))
  File "/home/genetic/.local/lib/python2.7/site-packages/CIRIquant/utils.py", line 82, in check_config
    globals()[i.upper()] = check_file(config['tools'][i])
  File "/home/genetic/.local/lib/python2.7/site-packages/CIRIquant/utils.py", line 49, in check_file
    raise ConfigError('File: {}, not found'.format(file_name))
CIRIquant.utils.ConfigError: File: /home/genetic/ncbi/dbGaP_26368/bwa_sam_CIRI/bwa, not found

###########
My "CIRIquant-Config.yml" content is as below:
tools:
  bwa: /home/genetic/ncbi/dbGaP_26368/bwa_sam_CIRI/bwa
  hisat2: /usr/bin/hisat2-2.2.1
  stringtie: /usr/bin/stringtie
  samtools: /usr/bin/samtools-1.9

reference:
  fasta: /home/genetic/ncbi/dbGaP-26368/fasta_file/hg19.fa
  gtf: /home/genetic/ncbi/dbGaP-26368/gtf_file/gencode.v19.annotation.gtf
  bwa_index: /home/genetic/ncbi/dbGaP-26368/bwa-sam-CIRI/bwa-0.7.17/hg19
  hisat_index: /usr/bin/hisat2-2.2.1/hg19

I apprecite it if anyone share his/her comments with me.

Do you see bwa under /home/genetic/ncbi/dbGaP_26368/bwa_sam_CIRI?

Dear @haci , I installed bwa aligner in bwa folder. In fact, in "/home/genetic/ncbi/dbGaP_26368/bwa_sam_CIRI/bwa" bwa is a folder and the executive file of bwa is in  bwa folder.  should I Change this path structure?

I believe adding /bwa right after bwa in your tools: section of the yaml file is worth a try.

Dear @haci, I did your recommendation in tools :  section and was passed but Now I have the same error for refrence :  section for fasta file as below:        CIRIquant.utils.ConfigError: File:/home/genetic/ncbi/dbGaP-26368/fasta_file/hg19.fa , not found . I appreciate it if you share your comments with me

@haci I downloaded hg19.fa from : [link] (https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/)     through this link:   [link] (https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz)

Once again, do you see hg19.fa under /home/genetic/ncbi/dbGaP-26368/fasta_file/? h19.fa.gz is not the same as hg19.fa, the former is compressed. Most aligners if not all are capable of working on compressed fasta files so you can give a try by altering reference section in your yaml file to reflect that.


text:
semantic similarity measurement for cell line ontologies
I have a set of cell line pairs and I want to know to what extent the pairs are similar based on their ontologies.
The problem I have is that I have found a Python library called Fastsemsim, but it has many similarity measurements. Or to be exact, there are actually many ontology similarity measurements. And I want to know which similarity measurement is best suited for my case (cell line ontologies similarities). There might be a couple of good ones that are good for cells ontologies.
These are a set of similarity measurements that are available in Fastsemsim.
Semantic Similarity

SemSim.TermSemSim
SemSim.ObjSemSim
SemSim.ObjSetSemSim
SemSim.SetSemSim
Specific Semantic Similarity measures
SemSim.MixSemSim
SemSim.avgSemSim
SemSim.maxSemSim
SemSim.BMASemSim
SemSim.SemSimUtils

Source:
https://fastsemsim.readthedocs.io/en/latest/library.html

Welcome to the site. Which is your question. What is the best semantic similarity? Or something else (do you want to compare the cell lines with non-semantic similarities measures)? Do you have other information besides cell lines ontologies (not sure were do you take those from or in which format though)?

The format is OBO but that is not related to the question. My question is which similarity measure (listed) is best or does better at finding semantic similarities between cell lines

About the metric you use, there are several reviews on the topic, if I recall well I think that the information content based are the best ones, of those, I only recognized the JiangConrathSemSim listed on the webpage.
However, if you want to use a semantic similarity measure I would use the BMA method, which is more comprehensive than other methods. But the BMA method is only a way to aggregate several metrics of ontologies into a single numeric value. I am not confident that FastSemSim implements correctly those methods.
Note, that according to the license, you should contact the authors before using it:

You must contact us before using it or its results or any application/tool based on it in any published work.

I am not sure that the cell lines fit in a ontolgogy. Because ontologies are defined as directed acyclic graph (DAG) and there are some cycles that happens in cells like the induced pluripotent stem cells show. However they might be right if you just compare the cell lines and not the state of the cells on those lines (use it only if you want to compare the lineage of your cells lines).


text:
What is the perferred method of optimization or energy minimization of small molecules downloaded from PubChem?
I will be docking a small library of molecules that I have downloaded from PubChem using AutoDock. I thought why not minimize their energy before docking using ChemDraw. I am thinking; what is the best method of optimization/ energy minimization of small molecules. MM2, MMF94 or other? And why?

The difference between different force-fields is not going to be major, it is the side steps which are.
When
If you are starting from a SMILES string, optimisation is a must obviously.
If you are using a 3D conformer from PubChem or even an actual sub-1 Å crystal structure from CSD, optimisation is nice for consistency.
Which
MMF94 is a solid choice. RDkit offers MMF94 (Merck force field 1994) and UFF (universal force field). The latter has some extra features, such as not kicking the bucket at atoms that aren't your usual suspects. ChemDraw offers MM2 and MMF94 IIRC, the former is very old. g16 in Gaussian (paid software) uses Gaussian Electrostatic Model, which is polarisable and will solve partial charges in a nice way (which I'll get to in a minute). OpenBabel uses MMF94 under the hood. You can even optimise with macromolecular forcefields, such as AMBER and CHARMM.
Given that the AutoDock uses a flexible ligand it should not matter too much anyway.
What else?
Autodock does not require conformers by virtue of the flexible ligand business (the sampler does rototorsions). Other programs might if they use rigid ligands (e.g. Fred) or they use a mixed approach and use a conformer library for speed (e.g. Rosetta ligand_dock). Conformer diversity is actually a bigger issue as conformers that are only slightly contorted, often even less than kB (0.6 kcal/mol), may not be sampled.
The thing that messes things is incorrect protonation. OpenBabel can protonate at a given pH, RDKit cannot. Not an issue for ligands with drug screens (cf. Lipinski rule), but is for natural substrates. However, for the protein it is essential to make sure acidic and basic residue are correctly protonated and that histidine have the correct tautomer —there are tools, c.f. MOE, that can do it but by far the best approach is reading up on the mechanism.
The next thing is dodgy Gasteiger Marsili partial charges —for Autodock 4 (and most other full/hybrid physics-based docking programs), partial charges are important: after all PDBQT has a Q in it for charge! For prepare_ligand4.py you give it a mol2 file. A major consideration is what format you are using for your molecules. This sounds a very silly thing to point out, but you'd be surprised how many trip up. An mol file (or an sdf file, which is the same but multi-entry) does not have an official partial charge column (although they are often added as comments after M END). So whatever force-field you chose make sure that the program you use gave you partial charges are not formal charges. This does not quite apply to AutoDock Vina, which is charge independent and more reliant on machine learning derived factors.
Verdict
Open babel is nice in that it can start with a SMILES, handles multiple formats including mol2 and PDBQT, can protonate at pH 7, uses MMF94 and does Gasteiger Marsili charges (and its confab algorithm is very good, not need here though).


text:
Do I need to study alternative splicing and isoformswitch separately?
I have done a study on isoform switch between different tissue type. Now, many people asking me to study alternative splicing also. But, dont they will give me same results? If there is a significant splicing difference; that will definitely leads to isoform switch. Or my understanding is incorrect ?
Thanks

Imaging isoform A uses exon 1, and isoform B uses exon 2. Now if you do an isoform switch analysis, you might see that in one condition isoform A (and exon 1) is preferred over isoform B (and exon 2). You might get similar results for alternative splicing. But if there is another isoform C that uses both exon 1 and 2, then the alternative splicing analysis of these exons are going to be affected by the expression level of isoform C as well.

Thanks @PhoenixMu


text:
Calling isoforms from long read data generated from partially degraded RNA
What will be the best tool to call isoforms from long read data generated from partially degraded RNA. By mistake we processed some samples with poor quality RNA to generate long read. Now we are looking for some tool which can merged the reads for the potential degraded transcripts. I think majority of the available tools like Squanti, Talon, bamboo, flair considered independent read as full length. Do the community has any suggestions?
Thanks

I haven't used these long-read RNA-seq tools so I could be wrong here. Transcript reconstruction combines reads from the same transcript and should be robust to partial transcripts. In addition, Nanopore direct RNA-seq produces a lot of partial transcripts. I tend to believe many transcript reconstruction tools can deal with them.
As to what tool to use, StringTie is a popular one and has been shown to be better than Flair. IsoQuant seems promising. Its developers argued it is more accurate than StringTie, Bambu, Talon and Flair. Squanti only works for iso-seq.

You could have a go using LAST. I've been developing and using a protocol for mapping cDNA reads that is based around LAST-mapping since 2017. The most recent RY model addition in LAST seems to clean up the majority of the random gene counting noise issues that we were having.
See the gene counting part of the protocol here. If you're working from direct/native RNA rather than cDNA (or don't care about strand-specific mapping), then the strand correction in previous steps is not needed:
https://dx.doi.org/10.17504/protocols.io.5qpvonn2bl4o/v17


text:
Single Cell RNA seq Analysis for low input Cells?
I am having a single cell RNA seq data from ~500 Cells. I do understand this number is very low for present day analysis tools like Seurat. So, I want to know what if one have a such a small number of cells which tool works best?
Thanks

you can perfectly use seurat as part of your analysis pipeline. It is hard to suggest tools or statistics without knowing the type of data you are looking at and what questions you want to address.

@Mack123456 Thanks for your comment. I am trying to use single cell RNA seq data on MiNION. On a trial basis I just started with ~500 Cells.

Hi, you will need to forward as much information as possible, or simply accept @Mark123456 answer. The problem is that we don't know the objective of the study.

@MichaelG. Thanks. I am ready to accept the comment as answer. But I cant find the accept answer tag in comment section.

Hi @user3377241 no need these are just comments not a full answer. If you look at the other questions it's detailed technical advice rather than broad answers. If you have a specific technical issue you can return and ask about it

You can perfectly use any of the single cell processing packages out there. Seurat, Monocle3, Scanpy, Kallisto, etc. Following some of the tutorials and vignettes will most likely help you figuring out how to proceed.
These packages provide useful functions to explore single cell data, process the data and each of them has a series of advanced features to address specific questions. You don't need 10k cells to look at number of gene/cell, transcript/cells, % mitochondrial genes, etc. You can perfectly normalize and standardize 500 cells and run PCA, and analyze the complexity that describes your data.
To address if 500 cells is enough for you pilot we need to know what question you want to address.
The following tool for example gives some intuition on how many cells to sample if you wish to detect n-cells per cluster.
https://satijalab.org/howmanycells
Hopes this helps


text:
"variable.features.n" in SCTransform
Whats the role of "variable.features.n" in SCTransform function of Seurat? I try setting "variable.features.n" at three different numbers (2000, 5000, 10000). In the downstream cells are clustering in different groups?

variable.features.n sets the number of features (you can think of genes in the case of scRNA-seq) you would like to use for the downstream steps such as clustering. Basically, using a gene that is expressed more or less at similar levels across different cell types would not be informative in terms of differentiating (for example via clustering) these cell types, this parameter helps you to choose the "most informative" (most variable across the whole data set) features.


text:
what does a sequence look like before alignment?
I am confused with Sequence alignment, which is

a way of arranging the sequences of DNA, RNA, or protein to
  identify regions of similarity that may be a consequence of
  functional, structural, or evolutionary relationships between the
  sequences.

to understand this easily, I would like to use this image (image_1)

the image above (image_1) is the primary part of nucleic acid structure (primary, secondary, tertiary, and quaternary), shown below (image_2)

is the sequence gcgc... an aligned sequence? if yes, what does the sequence look like before alignment?

Sequence alignment typically refers to the primary sequence. That is, we ignore the complex 3D structure of the molecule and its spatial arrangement in the cell. Instead, we use a symbolic representation composed of symbols like A, C, G, T, and U. Alignment is a conceptual tool that, as you pointed out, is used to find regions of similarity between sequences. Therefore, it only makes sense in the context of two or more similar sequences.
By way of example, consider a simple case of two sequences X and Y. This is what they would look like before alignment. The similarities between the two sequences are unclear.
X: ATTGACATAGGCACTTCGACTA
Y: GACATCTTCGGCTAA

This is what they would look like after alignment. The similarities (and differences) are now much easier to see.
X: ATTGACATAGGCACTTCCACTA-
      |||||     ||||  |||
Y: ---GACAT-----CTTCGGCTAA

This arrangement identifies regions of similarity that may be related through evolution. This doesn't always mean perfect agreement though. As this example shows, sometimes mismatches or gaps are allowed so that long stretches of identical matches are lined up correctly. Sometimes, for very closely related sequences, the optimal alignment doesn't require any gaps at all. This seems to be the case with some of your recent questions about sequence logos.
Note: Sequence alignment is a fundamental concept in bioinformatics, and there is a tremendous amount of research and development on the topic. It's hard to imagine doing any amount of serious DNA sequence analysis responsibly without a basic understanding of sequence alignment. I would recommend searching Google and/or Google Scholar for intro-level reviews of the topic.

There are a lot of really good materials by Ben Langmead (the guy behind Bowtie  and HISAT) on sequence alignment, what it means, and how it's done. I'd also recommend most things written by Heng Li, if you're looking for places to start learning about sequence alignment


text:
How to check if a given gene is expressed in a group of microarray samples if I do not have control group to compare with?
I have a microarray gene expression dataset consisting of placenta samples. I want to check whether these placenta samples are mixed with maternal decidua tissue. I have marker genes for decidua samples SCARA5, PIGF, SLPI. How to check whether these genes are expressed in placenta samples? In other words, how to check if a gene is up or down regulated without control group?

Which microarray? Some have negative controls that can be used to assess the probability that a gene/transcript is expressed.

I look at theese datasets: GSE73374, GSE73685, GSE9984, GSE122214, GSE22490, GSE35574, GSE37901. The microarray platforms are Affymetix and Illumina: affy_hugene_2_0_st_v1, affy_hugene_1_0_st_v1, affy_hg_u133_plus_2, illumina_humanwg_6_v2. Unfortunately I do not have access to illumina raw data, only to the preprocessed one. Can I still get negative controls from it?

At least some of those have the .CEL files included, so you have all of the data.

@DevonRyan I've looked .CELL files' description and indeed there are control probes. Could you give me a hint, how do I move from these raw data to the lists of up- and down-regulated genes?

Finding up and down regulated genes requires using R packages like affy and limma. Please search for tutorials.

Thank you @DevonRyan I will look them up more closely. But aren't these methods require having case and control groups of samples and not just one group like in my case?

Your previous comment asked for up and down-regulated genes. That requires some sort of treatment groups, yes.

Oh, I see now. I said up- and down-regulated, but I meant expressed and unexpressed relative to control probes. My bad.

You ask about which "genes are expressed" and then you mention "if a gene is up or down regulated". These are different, and given your application I think what you actually want to know if these marker genes are expressed. This is not a question of up/down-regulation relative to a control, so you do not need control samples for this.
It seems most of your datasets are Affymetrix, and there are R/Bioconductor packages that will give you the probability a gene is expressed. One such I found is frma (https://bioconductor.org/packages/release/bioc/html/frma.html), which returns an object "weight" that is "a vector of weights which roughly correspond to the probability of expression for each gene." I haven't worked with microarrays in a while, but I recall many other R/Bioconductor packages offering this functionality.
A detection p-value is calculated for Illumina array probes in Illumina's native BeadStudio software.


text:
How to get a list of all human morphogens?
Gene Ontology provides many genes annotated as taking part in various morphogenesis processes but I want to get a list of all morphogen coding genes specifically. Uniprot  does not have "morphogen" keyword.
I need some gene annotation database that has genes labeled as "morphogen".


text:
What is a PacBio "movie file"?
I came across references to "movie files" from PacBio sequencing in this paper:
https://www.jimmunol.org/content/204/12/3434
Specifically:

Movie files used to generate results presented in this article have been submitted to the National Center for Biotechnology Information BioProject (https://www.ncbi.nlm.nih.gov/bioproject/) under accession number PRJNA389440.

Looking at the SRA entries, it looks like it's just the raw reads as I expected.  But what is meant by "movie files"? (This is probably just my ignorance about PacBio as opposed to Illumina or others but I've never seen that term before and can't find much online.)

[disclaimer: not a PacBio person, I'm just vaguely aware of the technology]
Movie files are the PacBio raw data format, representing observed signal intensities in each sequencing well over the course of the sequencing run. The movie files can be used to re-call sequences with updated statistical models, which can improve accuracy, or test out hypotheses about how epigenetics alters the rate of nucleotide incorporation (i.e. methylation).
Illumina has a similar static version of this which is an image (or multiple images, depending on chemistry) of the flow cell after each synthesis cycle. Illumina's raw data format consumes tens of terabytes per run, and is of minimal value after base calling, so is almost always discarded after base calling.
Anything uploaded to NCBI SRA must contain FASTQ information. While the original uploaded files can be kept (which may include time/signal information), the files are converted into a FASTQ-like format for NCBI's internal representation. As such, there's a chance that movie files could be uploaded to SRA, but those files would also need to include the called FASTQ information (e.g. in a HDF5 container). In contrast to this, ENA does allow uploading raw data files without called information.

Got it!  In that paper's context do the authors probably mean they deposited reads in the SRA post-base-calling, then, rather than pre?  (I deposit FASTQ from my Illumina runs, not the raw images, for example.)  On the other hand I see NCBI's docs mention HDF5, so could it include signal intensity?


text:
When using cutadapt, can I specify an R2 adapter as optional when I have a required R1 adapter?
I'm using cutadapt 3.5 to trim adapters and perform some filtering on paired-end data.
Both R1 and R2 sequences have 3' adapters that might be found depending on the sequence length, but R1 also has an anchored 5' adapter that I want to require for the pair to be kept in the output.  For example, if that required R1 5' adapter is AAAAAAAA, the R1 3' adapter is TTTTTTTT, and the R2 3' adapter is GGGGGGGG, and I have this input:
R1:
@read1 should be kept
AAAAAAAACGTCCTGGGATTTGTAATAATATTTTAGTTCTGAGCGACAAGTAAGGGATAATTTTTTTT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read2 should be filtered
TACCAACTACATTTAGCTTCAGGCTAGTGATGCCCGCCGTCGGCACACTGGACACATGGCTTTTTTTT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3 should be kept
AAAAAAAACAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

and R2:
@read1
AGGCTCACCGCCACTGTTGTACCTTCCTATCGCCACTCTAAGATCTATGTAACCTCTCCCGGGGGGGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read2
TGAATAATACTGCGGTTTAATCGATATATTCGGGATTATGCAAGACACCCTACGTACTTAGGGGGGGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3
TACTGGAGGGTACAAGCCCGGACTATCCACGGCGTCAGCTGGCTTAGCATTGAAGGTCGGCGGTGTGT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

I'm aiming to get this R1 out:
@read1 should be kept
CGTCCTGGGATTTGTAATAATATTTTAGTTCTGAGCGACAAGTAAGGGATAA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3 should be kept
CAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

And this R2 out:
@read1
AGGCTCACCGCCACTGTTGTACCTTCCTATCGCCACTCTAAGATCTATGTAACCTCTCCC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3
TACTGGAGGGTACAAGCCCGGACTATCCACGGCGTCAGCTGGCTTAGCATTGAAGGTCGGCGGTGTGT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

This approach doesn't work as it discards read3:
cutadapt -a '^AAAAAAAA...TTTTTTTT' -A 'GGGGGGGG' example.R{1,2}.fastq --discard-untrimmed -o out.R1.fastq -p out.R2.fastq

I expect that's because the required/optional behavior only applies to the linked adapter for R1, and R2 is considered untrimmed for read3, which then filters the whole pair because of the default --pair-filter=any behavior (if either read is untrimmed, the pair is discarded).  I found I can get the filtering I want with --pair-filter first, but then the other filtering criteria like quality and length will only apply to R1, too.
My question is, can I get the trimmed/untrimmed behavior I'm looking for (read1 and read3 trimmed, read2 untrimmed) while also applying cutadapt's filtering to both reads in each pair? In other words, can I specify to cutadapt that R1's anchored adapter is required and the other two optional, without resorting to changing the pair filter setting?

This is not possible in a single cutadapt call, but can be done fairly seamlessly with a pipe between two calls and with cutadapt's interleaving feature.  (Thanks to the author for help with this.)

Giving the R2 3' adapter to one cutadapt command and interleaving the read pairs to piped output, a second cutadapt command can take that interleaved output and use both of R1's adapters as a linked adapter:
cutadapt -A GGGGGGGG --interleaved example.R{1,2}.fastq | \
cutadapt -a "^AAAAAAAA...TTTTTTTT" --interleaved --discard-untrimmed -o out.R1.fastq -p out.R2.fastq -

Now I can add other filters to the second call.
An additional read pair that should be filtered:
R1:
@read4 would be kept but R2 is too short
AAAAAAAACAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

R2:
@read4
CTACTAGACATATACCCTGGGGGGGGTTCCCCGTGAAGACGGTTACGATGTGGACAAGCACTCCAGCG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

I can now give a minimum length that applies to both in each pair:
cutadapt -A GGGGGGGG --interleaved example.R{1,2}.fastq | \
cutadapt -a "^AAAAAAAA...TTTTTTTT" --interleaved --discard-untrimmed --minimum-length 30 -o out.R1.fastq -p out.R2.fastq -

With that I get read1 and read3 in the output, as intended.


text:
Why do these NCBI representative genomes for ape species have no Y chromosome?
I recently downloaded a genome from NCBI for chimpanzee and was surprised to see no Y chromosome.  I added in one from a separate accession and carried on, but then today went to do the same for gorilla and saw the same thing, no Y.

For Gorilla gorilla (western gorilla): https://www.ncbi.nlm.nih.gov/genome/2156
For Pan troglodytes (chimpanzee): https://www.ncbi.nlm.nih.gov/genome/202

Why would this be?  I can see a separate entry for a Y chromosome for that same exact chimp used for the representative genome (Clint) here, but I don't see why it's not in the representative genome listing itself:
https://www.ncbi.nlm.nih.gov/nuccore/DP000054.3

tl;dr: technical difficulties, or sex
For the gorilla, that genome is from a female gorilla (Kamilah the gorilla) so she doesn't have a Y chromosome.
For the chimp link, there is a visualization of the chromosomes that you can click on that takes you to the genome browser for the Y chromosome (see image). Notably, there is a more recent genome build that does have a Y chromosome sequence assembly that is directly linked, they just haven't updated the main reference sequence yet.

To directly address the question, the Y chromosome is very hard to assemble. The easy part of it to assemble is already on the X chromosome anyways (the pseudoautosomal region), and the hard part of it is often a mess of repeats, especially in mammals. Many groups will decide to just assemble the X as well as they can and call it good, due to the difficulties associated with assembling the Y. If they are working with a male individual the Y sequence is still there, it just may not be anchored to a chromosome sequence, so it's only in the unscaffolded/unanchored assembly ("master WGS" I believe, usually). It is also probably highly fragmented (it isn't a whole chromosome of sequence stitched together, but rather a lot of little pieces of sequence).
In some cases they will focus all sequencing efforts on female individuals specifically to make their assembly job easier, because in females the X is at the same copy number as the autosomes and doesn't have interference from homologous regions on the Y.


text:
Use of Electronic Phenotype in EHR
May I know what's the use of Electronic Phenotyping using EHR data?
I did refer this link but have few questions
I understand that Phenotype is a set of criteria that you apply on your EHR data to select patients of interest from EHR database.
However what I am trying to understand is the difference between ordinary SQL query and electronic Phenotype
If both are same,
1) what I would like to know is the use of electronic phenotyping? How it is used?
2) Why do we have pheKB repository and other repositories? What's the use in validating them in multiple sites?
3) Is it like identifying patients real time? I am not able to get the usefulness of this. I understand that we select patients of our interest (cohort), but i can just write an ordinary SQL and send it to my EHR database to identify patients. but what's special about Electronic Phenotype?
4) I also see that they use multiple conditions (Lab tests, Medications, Diagnosis etc) to identify patients. What am I missing here?
Can someone help me by giving an simple example? I am new to healthcare and trying to learn. Would really be helpful

I understand that Phenotype is a set of criteria that you apply on your EHR data to select patients of interest from EHR database.

No, a phenotype is a way of behavior or other observed characteristic of a person, resulting from a combination on its underlying biology and the environment it's in. Think things like, height, running speed, strength, etc.
In the context in which you're interested, this would be more properly, "the clinical presentation of a disease or disorder".

Step one of diagnosing a patient is defining exactly how they present (i.e., exactly what their maladies, aberrant electrolyte levels, genetic mutations, etc. are).
An algorithm that only works on the test population is useless. Google "model overfitting".
The point of an algorithm is to take the EHR data and progress to a diagnosis. This isn't real time, you'd just run a batch of algorithms as test results come in.
You need patients in your test/validation sets to have the required inputs for your algorithm.

FYI, I don't know of any regulars on this site that deal with EHR data. If there's a community of such people on the web they must be elsewhere.

thanks for the response. Upvoted. Can you please help me understand this The point of an algorithm is to take the EHR data and progress to a diagnosis. ?

What is unclear about that?

I would suggest looking at the documentation for phenotips software and playing with the online demo.  This would show you some examples of phenotype data and how a database can be implemented to store it.  You may also want to familiarise yourself with HPO terms (Human Phenotype Ontology), a standardised ontology for recording human phenotypes (physical symptoms and characteristics which may be related to the underlying genotype).


text:
Appropriate value to replace outliers in clinical studies
I am working on clinical study where I am focusing on binary classification problem.
There are few variables which has outliers.
I am trying to know which would be the right value to replace outliers in a clinical study?
Is it mean, median min,max or mode of a column of interest?
Currently I use min and max for values below and above 1.5*IQR threshold
I use IQR approach to identify outliers but would like to know how it is done usually in clinical studies?
Can anyone help me with this?

In my experience, it is more common to use all clinical data as-is for clinical studies. And if data is missing, either omit the sample or omit the variable with missing data.
If your classifier can't handle the wide variation commonly seen in clinical studies then you may want to use a classifier which is less impacted by outliers.

Thanks for the response. Upvoted.

Removing outliers is common practice in statistical modeling and perfectly acceptable. However, with regards 1.5 IQR I am far from certain about this approach. Normally, if you want to be conservative then 3 standard deviations (SD) denote an outlier, which is more stringent than IQR. Some use 2 SD. If the value is lower than 2 SD from the group mean it isn't an outlier, that's because 1.96 SD is the normal distribution.
Normally you would delete the outlier and it takes no further part in the calculation, assigning it to 'na' should work. The subsequent calculation will recognise an 'na' value, I think in python it is np('na') ... Np = numpy.
In this case you appear to be deleting the value and treating it as missing data. Treatment of missing data is complicated I agree 100%. Mean, median are possibilities for missing data I can't see this working for binary data. One approach is regression and in python there is a command specifically to do this (which I forget), R will have the same command and it will assign a value  for the missing data. However, you can't do this for outliers.
The other routes to remove outliers are via transformation, log transformation is the most common. Transformation is a big subject.

thanks upvoted for the response


text:
Filtering raw sequencing reads
I have to filter the raw sequencing reads based on the following criteria:

Remove reads containing adapters
Remove reads containing N > 10% (N represents base that could not be determined)
Remove reads where the Qscore (Quality value) of over 50% bases of the read is <= 5

The reads are paired-end fastq files. I have been able to do the first part by using trimmomatic:
trimmomatic \
  PE in_1.fq.gz in_2.fq.gz \
  pair_out_1 unpair_out_1 pair_out_2 unpair_out_2 \
  ILLUMINACLIP:adapter.fa:2:30:5

The other 2 steps are a bit more challenging since I have not found a tool that perform exactly what is required. I could write some custom code of course but I would like to use the existing tools if possible. No need to reinvent the wheel.
Do you have any suggestions on how to proceed?


text:
removing nextera transposase adapters, cutadapt
I am learning about NGS analysis and im currently learning about QCing and removing adaptors.
I am working on SRR1972920_1.fastq file.
When running fastqc tool on that file, adapter contamination is present in the form of Nextera Transposase adapters.
Searching for that adapter sequence via google the sequence for at the 3' end to be removed is TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG (obtained via illumina adapter sequence-Nature).
Running a very basic cutadapt command line to remove the adapter
cutadapt -a TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG -o <any output name> <my input file>

After running the fastqc tool again on the new output file from cutadapt, it seems nothing changed and the contamination still present.
However, when running the same command above but changing the sequence to CTGTCTCTTATA, which is the sequence of transposase adaptors according to fastqc tool adapters_list.txt file all is good after running fastqc tool on the file
cutadapt -a CTGTCTCTTATA  -o <any output name> <my input file>

Now, why this is happening? which Nextera transposase sequence is the correct? or am I doing something wrong?

The adapter sequence you have googled is the sequence on the adapter primer. This works when you want to remove primer-dimers.
With transposase adapters or ATAC seq, you have very short fragments and you sequence into the adapter. So what you is sequenced is actually the reverse complement. See below image from trimmomatic manual

So if you look at what worked, and the sequence you googled, they are reverse complement of one another.
You can also see this in the NexteraPE-PE.fa, Trans1 and Trans1_rc @Haci recommended:
more NexteraPE-PE.fa
>PrefixNX/1
AGATGTGTATAAGAGACAG
>PrefixNX/2
AGATGTGTATAAGAGACAG
>Trans1
TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
>Trans1_rc
CTGTCTCTTATACACATCTGACGCTGCCGACGA
>Trans2
GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
>Trans2_rc
CTGTCTCTTATACACATCTCCGAGCCCACGAGAC

To answer your question, if most of your fragments are short, and there's not a lot of primer dimers, then yes you are doing the right thing.

Ok, im getting confused already :). What do you mean by adapter primer?

I meant the sequence of the adapter itself. Referring to the image I have in the post, adapter would be in the blue part. If you have short fragments that exceed sequencing length, you will sequence the adaptor on the right side(red). However the sequence you get in the fastq is the reverse complement.

I cannot quite help(*) on the problem you are having with cutadapt but can point you out to Trimmomatic, for which the developers have been granted permission to distribute  Illumina adapter sequences.
If you download Trimmomatic, you will see a bunch of Illumina adapter sequence files in the /trimmomatic-0.39/adapters/ folder. Actually the adapter sequence you have found via googling is listed in the file NexteraPE-PE.fa.
* Regarding the undesired output with cutadapt, you were probably removing one of the multiple Nextra adapter options but not the exact one used in your library preparation step. You won't have to worry about finding the exact adapter sequence used in the library prep step when using Trimmomatic as it scans for the presence of all possible adapter sequences.

Note that Trimmomatic has a trimming step that specifically looks for (and trims out) adapter read-through.

@Haci could you please elaborate how to let trimmomatic scan for the presence of all possible adapter sequences? I have just read the trimmomatic manual and I dont think I came across such option. Also I read on a blog post confirming what you saying "I know that a lot of people run Trimmomatic without checking which types of adaptors they actually have in the data and it often works. "

What I meant was the adapters and their reverse complements, PCR primers, ... are all included in the "adapter files" in the "adapters" folder. One should know at least the chemistry used as you have to specify an adapter file to the ILLUMINACLIP argument (not necessarily the exact sequences though). If I wouldn't know the exact chemistry, I would merge all the fasta files in the adapters dir into a single file and use this in the ILLUMINACLIP argument. This probably would not be necessary though, FastQC checks for/reports adapter sequences, that would give a hint on the chemistry used.

@Haci: that what I did initially yesterday but it did not work so I decided to write again here. Today I added the full path to the Nextera file and worked finally. My question is there a way to not write the full path for the file in the illuminaclip step? Also, I still see in the report loads of over represented sequences after the clip, what should someone do with those?

For your first question, you can place a copy of the adapter file in the directory from which you are running your command. For your second question, over-represented sequences are not necessarily arise from adapter sequences. You might want to blast them to see what they might be (low complexity regions, contamination, ...) But if you really would like to remove them, just add them to the end of your adapter file as a fasta entry and Trimmomatic will take care of them.

And somewhat unrelated but since you are dealing with a lot of QC, I strongly recommend MultiQC (https://multiqc.info). It does not do the QC itself but gathers results from QC tools like FastQC and generates a nice and concise QC report of your samples.

@Haci, your help is invaluable. Much appreciated. Last question regarding Trimmomatic. In the manual it talks about leading and trailing option by removing the base based on a quality value you decide. The examples are showing LEADING:3. What does the 3 value represent here? quality phred score of 3 is  really low? or am I missing something?

This is from their website Remove leading low quality or N bases (below quality 3) (LEADING:3). Yes, a score of 3 is very low but I don't know why they decided on this default value. Would be interesting to ask as a separate question!

I came across the same problem of Nextera transposase contamination in my shotgun metagenome sequence. I specified the library in trimmomatic and Nextera transposase adapters were successfully removed. Please see the code below:
java -jar /opt/software/Trimmomatic/0.39-Java-1.8/trimmomatic-0.39.jar PE -phred33 1004_R1.fastq.gz 1004_R2.fastq.gz ../QC_data/1004_R1.fastq.gz ../QC_data/1004.qcup_R1.fastq.gz ../QC_data/1004_R2.fastq.gz ../QC_data/1004.qcup_R2.fastq.gz ILLUMINACLIP:/opt/software/Trimmomatic/0.39-Java-1.8/adapters/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

The output showed completed successfully with the following message:
TrimmomaticPE: Started with arguments:
 -phred33 1004_R1.fastq.gz 1004_R2.fastq.gz ../QC_data/1004_R1.fastq.gz ../QC_data/1004.qcup_R1.fastq.gz ../QC_data/1004_R2.fastq.gz ../QC_data/1004.qcup_R2.fastq.gz ILLUMINACLIP:/opt/software/Trimmomatic/0.39-Java-1.8/adapters/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
Using PrefixPair: 'AGATGTGTATAAGAGACAG' and 'AGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTCCGAGCCCACGAGAC'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTGACGCTGCCGACGA'
ILLUMINACLIP: Using 1 prefix pairs, 4 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 18205918 Both Surviving: 10467310 (57.49%) Forward Only Surviving: 7470543 (41.03%) Reverse Only Surviving: 48503 (0.27%) Dropped: 219562 (1.21%)
TrimmomaticPE: Completed successfully

And the QC report showed the Nextera transposase was successfully removed.
[![QC report][2]][2]


text:
Fastq-dump script download X spots or all
Im trying to write a script where an optional input of -X flag can be used, or if that info is not available download all reads.
my script as follow:
# usage : bash dump_qc.sh <run id> <number of reads>

#run id variable.
SRR=$1

#maximum reads variable
X=$2

fastq-dump --split-3 -X ${X} ${SRR}

So, if I run the script
bash dump_qc.sh SRR1234 100

Will download 100 reads from the specified SRR file.
What I would like to do is to make the second argument optional, if not specified download the whole data set!

You don't need to copy your arguments into variables, you can just use $1 directly. Anyway, this should be enough:
#!/bin/bash
# usage : bash dump_qc.sh <run id> <number of reads>

if [ -z $1 ]; then
  echo "ERROR: at least one argument need to be given!" >&2
  exit 1
fi

if [ ! -z $2 ]; then
    x="-X '$2'"
fi

fastq-dump --split-3 "$x" "$1"

Thanks terdon for your solution. Im really beginner in bash programming. Could you please be more verbose and explain your code? :)


text:
Unable to read all alignments in stockholm format : UnicodeDecodeError: 'utf-8'
I have a file in Stockholm format and I am trying to read it using the Biopython package. The code is as follows :
import numpy as np
from Bio import AlignIO

alignments = AlignIO.parse("Rfam.seed", "stockholm")
for alignment in alignments:
    print(alignment)

It gets read for some alignments and then throws up following error :
Traceback (most recent call last):
  File "test.py", line 5, in <module>
    for alignment in alignments:
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Bio/AlignIO/__init__.py", line 394, in parse
    for a in i:
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Bio/AlignIO/StockholmIO.py", line 408, in __next__
    line = handle.readline()
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/codecs.py", line 322, in decode
    (result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xe1 in position 2468: invalid continuation byte

How to resolve this error?

Can you not just remove the unicode character from your file? The error is due to a Unicode character (å)...

The last line indicates that this is an encoding problem. Try to see if the function you are using to read the data has an "encoding" parameter that you can set to whatever encoding your file is using. If that doesn't work, you can change the encoding of your file by:

with save as using a text editor of your choice and selecting "UTF8" as the encoding. That would be a quick work-around if you have a small number of files.

or

You can use the codecs library to check and convert file encodings within your code. Here is an example.

Warning: I haven't used Biopython in ages but encountered similar problems with Python in the past.


text:
pymol script to only select 5 chains within a distance from a reference and save the selection
I am generating 30 asymmetric units of a protein structure. Now, I want to select different combinations of 5 units which are 30 Angstroms away from the reference and save it as a different file. Here is the command I could work out till now :
fetch 3q2v
symexp sym,3q2v,(3q2v),30
select chains, bychain all within 30 of 3q2v

How to loop in and select only 5 chains at a time obtained from the setup obtained by the above code, and save it in separate files?

Firstly, I would preface that the PyMOL command line in the GUI is great, but you can run PyMOL as a regular Python library (calling the appropriate python). This makes life so much simpler for more complex manipulations. You can do really complex manipulations with it. I cannot suggest it enough.
The API version of the commands works only marginally differently, but you can do API calls in the GUI cmd.) and PyMOLWiki documentation is only slightly weaker for them (help() fills the gaps).
But if you really want to use the GUI route...
Let's get the names of the models in the selection
onames = cmd.get_object_list('chains')
print onames

Now, you want to save subsets of 5 chains (technically, objects or models). That is, you do not care about the order they are in and you don't want the same one twice (actually selections with 4 or less), so you want all combinations. Make sure to do the following
    print len(onames)
And google that number as <n> choose 5. That is how many files you will get. I am pretty sure most of these sets will be meaningful and will be useless.
In a script a loop is easy, while in the gui you have to keep it simple or use a list coprehension hackishly.
 from itertools import combinations
 for s in combinations(onames, 5): cmd.save('-'.join(s)+'.pdb', ' or '.join(s))

Obviously, you could alternatively use an enumerate method around the iterator and get an index for the file. But do note f-strings will not work.
 for i, s in enumerate(combinations(onames, 5)): cmd.save(str(i)+'.pdb', ' or '.join(s))


text:
SRA Toolkit execution problem
I am trying to install SRA Toolkit on MacOS, following all the instructions from NCBI website. It is installed but fastq-dump cannot be executed because "cannot execute binary files". I tried "chmod +x", noah from home-brew to make binary executable on MacOS but fast-dump still cannot be executed. I attached the screenshot of my terminal window. Can anyone help?

Did you upgrade to Catalina (a bad thing due to no support for any 32bit at all)?

You can't run Linux binaries on OSX (or vice versa). Download Mac OSX binaries.


text:
Common QC parameters and threshold in human WGS pipeline
I am trying to summarize a list of QC-parameters and its threshold that should be applied for human WGS in general to avoid bad quality of data.
Which qc parameters for WGS are you using and why?
Thanks!

Hi and welcome to the site! Could you please [edit] your question and narrow it down a little. We don't do discussions here, we're strictly Questions and Answers, so the more specific the question, the better. What parameters are you interested in? GC content? Coverage? FASTQ quality? Insert size distribution? What will the sequences be used for (different applications may need different approaches). Please narrow the scope of your question down.

I would suggest you look at the website for the excellent MultiQC tool (https://multiqc.info/).  We use it in a clinical setting to collate all the QC statistics from our pipelines into a single report (It displays results from multiple QC programs such as FastQC, Picard, bcl2fastq, VerifyBAMid etc).  The homepage has a selection of example reports - clicking on "Whole Genome Seq" shows a report with the metrics most labs would expect from their WGS pipeline.
If you dig around in their GitHub repo (https://multiqc.info/) you will get a good idea of the default thresholds set for flagging results - there should be a config file somewhere specifying all these values.

For your negative controls, you might want to check out QC Fail (https://sequencing.qcfail.com/).  There you can see examples of how things look in QC when things have gone wrong, such as loss of quality, position-specific errors on flowcells and contamination with primer dimers.


text:
Make a fasta out of mapped reads without taking into account the reference
I've aligned reads onto an haploid reference genome. I'd like to have a consensus sequence of all my aligned reads that doesn't take into account the reference genome I used (i.e. I just want the nucleotides from my reads at each base (taking into account frequency and quality) but if there is no read a one position I rather have N/? than the reference nucleotide.
Is there a way to do this ? Thanks

vSNP outputs a VCF file that includes positions with no coverage, represented as "N".  Then there is a tool to merge a VCF file into the FASTA reference  vsnp_merge_vcf_into_fasta.py.  It will merge those "N" regions into the FASTA reference. There are options for considering frequency and quality.
 It's limited to just SNPs and those regions along the reference with no coverage.  It does not merge small indels or unmapped reads.  Pilon is a better tool for applying indels, but I'm not aware of it letting regions with no coverage be represented as "N".
vSNP_step1.py -r1 *_R1*fastq.gz -r2 *_R2*fastq.gz -r ref.fasta

Use *zc.vcf output to merge into FASTA
vsnp_merge_vcf_into_fasta.py -f ref.fasta -v *_zc.vcf

This is similar to this question, which links to this solution.
The difference here is that you want to go through and replace low coverage (maybe zero-coverage?) positions with Ns rather than using the reference. You can likely accomplish this using samtools depth to measure coverage and then writing a script that will go through each position in your reference, check its coverage, and string-replace with an N. Probably a bash one-liner could do this if someone very devoted to bash one-liners came up with one (I don't care for them, but other people are more creative than I). Very quick and dirty pseudocode might look like this:
### Someone could come up with something better if they spent >3min on it.
### should be easy to get it working if you use biopython
# read_and_preprocess() stands for whatever you need to do to get data in the right shape
coverage_data = read_and_preprocess('samtools_coverage.txt')
# reference_with_variants.fasta comes from the solution linked above
reference_fasta = read_and_preprocess('reference_with_variants.fasta')

threshold = 0
fasta_string = ''
for position in reference_fasta:
  if coverage_data[position] <= threshold:
    fasta_string += "N"
  else:
    fasta_string += reference_fasta[position]

# nonexistent fn that writes a fasta.
write_as_fasta(fasta_string)

I am a little skeptical of this approach because it seems like if you do not trust the reference sequence to stand in for zero-coverage areas, then you probably shouldn't trust it to align to and call variants against. Read mapping can be a rather noisy process, especially when you are using a more distantly related reference. I think this is related to the @swbarnes2 comment.
I would suggest instead indicating the regions of low/zero coverage by putting them in lowercase acgt with your covered regions uppercase ACGT, which preserves the coverage information that you care about but still allows you access to the reference information. It's not a perfect solution, but I think it makes the most out of the marginal data you're working with, and it is trivial to change the above example code to do it.
Hope that helps-

Have you googled de novo assembly?

I can’t do a de-novo assembly as I have really low coverage. I’m working with really short reads from ancient DNA.

You can make a consensus fasta from your bam, but that uses the reference.  If your read depth is that low, can you be confident of discrepancies between your reads and the reference?

No, but that’s another problem. Thanks. The thing I want to do can also be useful when you want to reassemble only some targeted genes from medium coverage data I guess.


text:
Difference between isoforms and paralogs in transcriptoms?
I've assembled RNAseq data into a transcriptoms using Trinity. There's a option to keep only the longest isoforms for each transcripts and it lead me to wonder how it deals with duplicated genes (-->transcripts). I haven't found an anwser in the original paper. Althought I reckon paralogs and isoforms are two distinct features I don't see a way to make a distinction when you only have trasncripts.

It depends. If you are performing Illumina based RNAseq, which is likely, it is definitely tricky, You should be able to do this in an ideal world because,

a paralogue will contain SNPs, i.e. there will be nucleotide variation between paralogues, but the mRNA can of an identical length
isoform will also have identical sequence but the length of the mRNA will vary depending on what initiation codon it used

I sort of agree however because,

there will be 'grey-zones' where paralogues are different sizes and VCF/ SNP calling is not perfect and ...
assembly of paralogues can be difficult, risking collapsing a paralogue into a single mRNA, but it is not impossible depending on the level of SNP variation
assembly of isoforms would be really hard

I personally wouldn't know how point 3 is solvable, but I guess it is if your isoform is relatively small, i.e. well within the read-length of Illumina output.

If a reference genome assembly exists, then reference-guided assembly with Trinity will keep paralogs, because it treats each genomic region as a different assembly set.
If a reference genome doesn't exist, paralogs will probably result in the construction of non-representative "isoforms" containing variants present in multiple genes. Inchworm does a greedy assembly of reads, so if a chimeric transcript can be constructed from the reads (which is likely with paralogous genes), then it might end up being fed to the subsequent steps. Whether that is actually emitted at the end by Butterfly will mostly depend on how the paths through the putative assemblies are traversed, which probably depends on the tolerance of Bowtie to errors.
If paralogous transcripts are sufficiently similar that the differences are within the error tolerance for assembly, then I would expect that they'd be collapsed into a single transcript.

Thanks @gringer its worth doing assemblies both ways, IMO it would work well if your 100% on all paralogue and isoform possibilities. In organisms with massive multigene families (where you don't know the full repertoire of sequences) using de novo assemblies is doable to prevent the collapse of paralogues depending on the distance between mutations and is particulalry true now Illumina now has pretty long read lengths. In the days of the 150 -200 bp reads then yes collapose would happen.  Clearly minority paralogues can be lost regardless.


text:
How to get coding sequences (CDS) of UniProt reference proteomes?
I downloaded the following reference proteomes:
https://ftp.ebi.ac.uk/pub/databases/reference_proteomes/QfO/
'Drosophila':'UP000000803_7227.fasta',
'Arabidopsis':'UP000006548_3702.fasta',
'C.elegans':'UP000001940_6239.fasta',
'Zebrafish':'UP000000437_7955.fasta',
'human':'UP000005640_9606.fasta',
'mouse':'UP000000589_10090.fasta',
'yeast':'UP000002311_559292.fasta'

Now, I would like to have the corresponding coding sequences (CDS). My temporary solution is to download CDS files from ENSEMBL FTP server, translate those and match with the ref proteome. However, this only covers 50-90% of protein sequences.
Thanks!

I went through diverse README files and found the answer:
Coding DNA Sequence FASTA files (*_DNA.fasta)
=============================================

These files contain the coding DNA sequences (CDS) for the protein sequences
where it was possible.
The format is as in the following example (UP000005640_9606_DNA.fasta):

>sp|A0A183|ENSP00000411070
ATGTCACAGCAGAAGCAGCAATCTTGGAAGCCTCCAAATGTTCCCAAATGCTCCCCTCCC
CAAAGATCAAACCCCTGCCTAGCTCCCTACTCGACTCCTTGTGGTGCTCCCCATTCAGAA
GGTTGTCATTCCAGTTCCCAAAGGCCTGAGGTTCAGAAGCCTAGGAGGGCTCGTCAAAAG
CTGCGCTGCCTAAGTAGGGGCACAACCTACCACTGCAAAGAGGAAGAGTGTGAAGGCGAC
TGA

So, EBI itself provides the corresponding coding sequences, which they name _DNA.fasta. I misinterpreted the file-name as genomic (unspliced) DNA.


text:
Why my bim file doesn't match to my ped file as the Plink documentation suggests?
Plink documentation about .bim says the columns 5 and 6 are the Allele 1 and Allele 2 respectively.
The documentation about the .ped in turn, says "...The first six fields are the same as those in a .fam file. The seventh and eighth fields are allele calls for the first variant in the .map file ('0' = no call); the 9th and 10th are allele calls for the second variant; and so on.."
I've simulated a dataset as per:
$ plink1.9 --simulate ds1.sim acgt --out ds1
...
$ head ds1.bim
1       null_0  0       1       C       T
1       null_1  0       2       A       C
1       null_2  0       3       C       A
1       null_3  0       4       T       A
1       null_4  0       5       G       C
1       null_5  0       6       T       C
1       null_6  0       7       T       C
1       null_7  0       8       G       C
1       null_8  0       9       C       G
1       null_9  0       10      T       A

Than I converted it to plain PED and MAP:
$ plink1.9 --bfile ds1 --recode --out ds1
...
$ head ds1.ped > ds1.head.ped
$ vim ds1.head.ped
per0 per0 0 0 2 2 T T A C A A T A C C T C T C C C....

Given the data above, why my .bim file shows "C and T" for the first variant while my .ped file shows "T and T"?
Apparently they don't match as per described in the documentation (or more likely, I'm misinterpreting it)...
Also, how do those 2 files related to each other? E.g.: How do I know that my Sample 1 has the SNP's null_3, null_4 or null_5, etc?
I need to clarify those doubt because I want to fit a LightGBM model with the Genomic Data (can be plink or vcf). Not sure yet how to do. I was trying to convert the plink files to a CSV-style somehow, to be able to use in my model. Not sure if this is the way to go though...
Thanks!

Responded in the plink2-users group (https://groups.google.com/forum/#!topic/plink2-users/zgJxdXxvdLo ).


text:
What is the best way to compare treated/control samples?
I have two populations of mice, say RED and GREEN; there are 5 pairs treatment/control in each.
The table is the following:
gene_name  p-value(treatVSctrl) RED1_ctrl RED1_treat RED2_ctrl RED2_treat ... GRN1_ctrl GRN1_treat etc...
gene1      1.3                  4         10         5         10             6.2       9.75
gene2      2.8                  3         6          6         12             3         6
gene3      4.4                  5         12         4         8              7         10

I need to compare two groups to have an idea of the expression of different genes after treatment.
I also need to find a difference between RED and GREEN group in terms of genes expression
The thing is I cannot compare mice one by one just measuring the fold change, I need to consider p-value. What is the most correct way to perform this analysis?

Its an easy comparison but I was confused by the header p-values. I can't help there needs to be raw quantification. I am not aware of a stat that can take an unknown 'p-value' and compare to other p-values at a given level of significance

Thank you. Do you mean there must be a p-value for each pair?

No the header in the first column says " p-value(treatVSctrl..." . I don't understand the quantification, but the column header is stating this is a p-value.

Agree the "consider the p-value" instruction is confusing. The first column clearly doesn't represent p-values, as p-value are probabilities and must be in the range [0, 1].

@NuclearHoagie it could represent a percentage, i.e. (0 to 1.0) * 100 , e.g. this could be a pairwise T-test. Alot of values are hovering around the 5% mark. Where there are mixed results it would be a reason for the OP to ask how to move forward.

I don't really understand the table, particularly the 'p-value' in the header. More generally I don't really understand the experiment nor the data.
As a general rule the Wilcoxon would be a reasonable guess of the test using raw data. This is non-parametric version of the T-test without the rigid assumptions of parametrics. Control mice go in one column and subjects in the other. Scipy.stats has the test, its certainly in R (which I don't really use).
print (scipy.stats.wilcoxon(x = [1,2,3,4,5], y=[1,3,5,6,7], zero_method='wilcox', correction=False, alternative='two-sided', mode='auto', *, axis=0, nan_policy='propagate'))

There you go, filled out with some example data.


text:
How can I separate upregulated genes from downregulated on a boxplot?
I have a table of fold change shifts for two different groups. Some of them are upregulated, others are downregulated. I wonder if it is ok to subset only positive values and to boxplot them to show only upregulated genes and vice versa?
     Gene   REDmean   GRNmean
1    gen1   4.918239  4.547266
2    gen2   3.635476  3.447484
3    gen3   -3.759263 -3.283889
4    gen4   3.962326  2.818917
5    gen5   -3.474217 -3.123105
6    gen6   3.653160  2.907355

Try again, use separate boxplot columns i.e. upregulated vs. downregulatated that are created using a dataframe e.g. pandas. For pandas use a filter to create a new column of 'UpDown', i.e. a boolean of upregulated or downregulated, and then use the groupby('UpDown') command. Send that into matplotlib boxplot function and it should work. The 'Gene' column needs removing (I think). My pandas is too rusty to knock out the code. Someone here might show a solution in R, but if you attempt a pandas approach and post it, someone will probably run with it.

Basically you'll reorder and pivot the data using a dataframe using the up- down-regulated filter (I'm not sure if 'pivot' is a correct term). What I don't like is your use of the word 'mean', for the real thing you need the raw data ... but this would be cool for creating the filter. "'Gene' column needs removing" well thats a bad idea (thats your results), but for the plot it will cause it to fall over if its dragged in on a dataframe - thats what I meant in the comment above.

Thank you. As I have asked in my previous posts I am still thinking if this is a valid approach to compare the mean fold change between two groups. And once again, I wonder if I can compare their fold changes (FC), say, for gene1 in RED group I have 5 for control and 25 for treated sample and in GRN group (for the same gene1) I do have 2 for control and 10 for treated. For REDs I do have a 5X FC increase (from 5 to 25) and for GRN it is also 5X (from 2 to 10). For me there is no difference between RED and GRN groups based on gene1.

The idea of the 'REDmean' and 'GRNmean' is that I have multiple samples, so I have taken a mean control vs treatment ratio (that are presented in a table above). The question is if I have the right to do that, and if yes - can I distinguish them just by separating negative FC from positive?

AS A PLOT - I have no problem. As a statistical test ... I think we covered the stats and I personally would not favour change-fold ... but if is an approach to standardisation (its a separate question) that is different. However, others might agree with you, i.e. a test on change-fold. The operation you are describing here to manipulate and plot is a serious question that requires formal code in R or Python. "Try again" ... sorry this was a personal comment, my first comment got deleted so I "tried again", this could be misinterpreted however.

It is good to keep questions distinct, because this question over a very large dataset isn't trivial, i.e. to provide the code that will perform the pivot. I honestly can't remember fully how to do it without looking it up (I've done similar things in the past).

I'll be working on pandas again shortly (couple of weeks). I'll be able to answer the question then if someone doesn't do this before.

Great! Thanks a lot for the help

Let us continue this discussion in chat.


text:
How can I get a column with real gene names in my ballgown analysis?
I am doing RNAseq analysis and I am using ballgown procedure for it.
There is an option in R to calculate differentially expressed genes and use FPKM in calculating differential gene expression, which is the following:
results_genes = stattest(bg_filt, feature="gene", covariate="state", getFC=TRUE, meas="FPKM")

When accessing head(results_genes) there is a column named id that contain gene names.
While some genes contain more or less readable names, others look like gene-LOC100037748, gene-LOC100037751, gene-LOC100037760 etc...
Can I create a column with real gene names that correspond to these? Is there a sort of API or something maybe?


text:
How to input the opposite strand of DNA?
I am trying to write a program that will count the repeats of a given DNA strand:
I am using a Counter library, something like this code shown below
from collections import Counter

        s="XXXXXAGCCTGCCAAGCAAACTTCACTGGAGXXXXXTGTGCGTAGCATGCTAGTAACTGCAXXXXXTCTGAATCTTTCAGCTGCTTGTXXXXXTGGGCCTCTCACAAGGCAGAGTGTCTTCATGGGACTTTGATATTTATTTTTGTACAACCTAAGAGGAACAAATCCTTTGACACTGACAAATTGGCTTCCATATTTTATACCTTAATCATCTCCATGTTGAATTCATTGATCAACAGTTTAAGAAAAAAAGATGTAAAAATGCTTTTAGAAAGAGAGGCAAAGTTATGCACAATAACTTCTCATGAAGTCACAGTTTGTTAAAAGTTGCCTTAGTTCACAATAAATAATTATGTATGCTCTATAATTTCAGTGA"

def chunk_string(s, n):
    return [s[i:i+n] for i in range(len(s)-n+1)]

counter = Counter(chunk_string(s, 5))

Repeat supposed to be "declared" as an identical sequence on the same strand (like XXXXX here that occurs 4 times). However, I can also have this repeat (for ex. XXXXX in this particular test case) on the OPPOSITE strand.
How am I suppose to intake the opposite strand in case of using .fasta files? As I know for today - fasta is presented only in a single strand way. Or am I completely mixed the terms?

Yes, you're going to want to take the reverse complement of your string in order to analyze the other strand, that is the standard.
Also, it is worth pointing out that you may want to use a generator expression instead of a list expression for your chunk_string function:
def chunk_string(s, n):
    return (s[i:i+n] for i in range(len(s)-n+1))

This will result in a generator that yields kmers as the Counter constructor iterates through them, which means that you don't need to first construct and store the entire list of kmers first and then pass to Counter but instead can feed them directly to it. Not essential for small strings, of course, but a good practice and useful for larger strings. See this SO post for a better explanation.
You can also use Biopython's reverse_complement function in order to easily obtain the reverse complement.

What's going to be simplest is to generate the reverse complement of your short search sequence, and search for that too.

Hi @swbarnes2, this looks like the start of a good answer; any chance you could expand on this approach?


text:
How to build advanced boxplots from datasets which columns contain only numbers?
I have got a fold-change dataset and I simply would like to plot something like this:

The problem is that I have only numbers in my columns and when I use a boxplot() function I can plot only boxes but not dots, jitter, etc.
I use this code to build a df and then proceed to a boxplot
A <- c(-1, 25, 3, 4, 50)
B <- c(23, -41, 32, 68, 26)
C <- c(63, 71, 22, -8, 2)
D <- c(16, -49, 2, 56, 6)
df <- data.frame(A, B, C, D)
print (df)
boxplot(df)

How can I plot something that looks like the above picture? Thanks

You want ggplot2, specifically the geom_boxplot and geom_dotplot functions. Check out the cheat sheet (warning: PDF download) to see the full functionality of the package.

Thank you. The problem is that I don't have geom data - if you check the dataset I have proposed - there is no "y" data, I have only 4 columns (A,B,C,D) and rows with numbers. It is ok for default boxplot() while I cannot add geom, can you specify what should the ggplot2 script look like for the above-mentioned dataset? Thanks

To coerce your data into a form interpretable by ggplot2, you can use the gather function from tidyr.
gather has been superseded by pivot_longer.
library(tidyverse)
df <- tibble(A = c(-1, 25, 3, 4, 50),
             B = c(23, -41, 32, 68, 26),
             C = c(63, 71, 22, -8, 2),
             D = c(16, -49, 2, 56, 6))

tidyr::pivot_longer(data = df, cols = c(A,B,C,D)) |>
  ggplot(mapping = aes(x = name, y = value)) +
  geom_boxplot() +
  geom_dotplot(binaxis = "y", stackdir = "center") +
  theme_classic()


text:
Is there a software that allows to convert a jpeg graph into dataset?
I have a very old graph as a .jpeg picture that I would like to convert into a dataset to be able to operate with it using R/python.
I have no idea if it is even possible.

How would that even be possible? How many datapoints are actually there? In theory, you could extract an infinite number, but the original data only had values for a small subset of all possible points on that line. Given that there are no visible points in your graph and that the quality is extremely bad, any extrapolated dataset would be a wild guess at best.

Looks like Graph Grabber can do the job.
https://www.quintessa.org/software/downloads-and-demos/graph-grabber-2.0.2#download


text:
How to compare differential expression of same genes between two groups?
I need a statistically significant "decision" that DEGs from group A are different (or not) from group B (after drug treatment).
For that, I have a mean fold change delta, calculated as the mean of 5 different treated samples minus the mean of the same, but untreated samples.
        A.delta    B.delta
gen1    4.9182385  4.5472662
gen2    3.6354761  3.4474844
gen3    3.7592630  3.2838886
gen4    3.9623264  2.8189166
gen5    3.4742171  3.1231050

In my vision, there is not much difference in fold change delta between A and B for the gen1 for example (4.91 is not very far from 4.54). Same for gen2. While for gen4 there might be something interesting.
What kind of plot may be the most useful for that kind of comparison?

I do feel we've been here twice before with this question https://bioinformatics.stackexchange.com/questions/18732/what-is-the-best-way-to-compare-differences-in-deg-expression-between-two-popula and this one https://bioinformatics.stackexchange.com/questions/18660/what-is-the-best-way-to-compare-treated-control-samples both appear very similar questions concerning mouse models. Originally labelled RED and GRN

If the answers you got before not satisfy you, you can always leave there a comment, or edit your question and clearly indicate what about the problem is that you can't use the proposed solutions.

This data is genuinely paired data and you have stated values are a mean. Thus the raw data can be analysed in your previous questions on the same theme described here:
What is the best way to compare treated/control samples?
And here:
What is the best way to compare differences in DEG expression between two populations?
If this is 'ANOVAerable' (its not a word, but if ANOVA is appropriate) than a two-way ANOVA would provide all the answers (however see below regarding the trend in the data). This analysis will perform a 'A' versus 'B' and each 'gen' compared to the other, as well as the interaction between A/B and 'gen'. You would need the raw data, you can't do anything with means (unless you are assuming the distribution). The distribution within each 'gen' is important (if its normal - conforms to the normal distribution, you're in business).
There is nothing further to add and I do feel we need to draw this series of questions to a close.
A regression model could be explored because there is a clear trend in both data sets, but I'm uncertain it will be significant. Within ANOVA however this might trap the outlier, but that is more complicated (you standardise against the trend thus 'gen1' is comparable to 'gen5' and therein 'gen4' can be identified as an outlier). This was touched on in the answer to the previous question (ols and python)


text:
How can I find all the columns available in UniProt.ws package for R?
I am trying to find the subcellular localization of my 10,000 proteins using UniProt.ws package for R. However, I am unable to find all the columns available for query.
I used another package named Bioconductor which had easy way to lookup all the columns available using COLUMNS(org.Hs.eg.db)
I tried reading the documentation and couldn't find anything which can output subcellular localization. I am new to this field, any help is appreciated.

From the UniProt.ws manual:

columns shows which kinds of data can be returned for the
  UniProt.ws object.

So, you will need to go over the output of columns(your_uniprotws_object) if it includes the data type you are looking for.
And here is a related question.

Thank you. It worked for me now with this code:
up <- UniProt.ws(taxId=9606)
keys <- c("4957","81551")
columns <- c("PDB","GENES","SUBCELLULAR-LOCATIONS")
select(up, keys, columns, kt)

Do you have a small sample of the proteins you are looking for and what you have tried so far in R using Uniprot.ws ? (including output of it), it could be helpful to identify your problem and how to solve it.
Personally, I would rather consider Bioconductor as a platform for providing R packages dedicated to bioinformatics than a package to analyze data (but it's my opinion, maybe I'm wrong).

I have 10000 proteins which I am looking for sub-cellular localization. I tried to find the columns available in UniProt.ws using COLUMNS(UniProt.ws) but it didn't work. Finally, I looked into Bioconductor as you said and found the answer. https://support.bioconductor.org/p/78380/ . res <- select(up, "P35579", "SUBCELLULAR-LOCATIONS", "UNIPROTKB")


text:
Should I use log2-CPM values (voom-limma) as input for my model?
We have created a model to integrate several OMICs data, but we realized that the maximum TPM values of RNA-Seq data were so big that had unexpected effects on our results. We hypothesized that this might be due to the negative binomial distribution inherent to this OMIC's layer.
To tackle this circumstance, we proposed to introduce as input the values obtained after using voom() (paper) in our count data. The reason behind that decision was that voom would transform the input data to be more accurately introduced to a linear model. There are, though, several issues/arguments against doing so:

The results obtained through voom() contain negative values, which do not make sense in our model. I tried to correct that by adding to all the results from a specific sample the absolute minimum value of that sample, but I am not confident whether that is a good way forward.
Original 0 values (non-expressed) turn into different negative values across samples (-6 in Sample 1, -4 in Sample 2), creating a difference that did not exist to begin with.
Some sources do not recommend to log-transform count data [1], including data from negative binomial distributions [2].

I am not really versed into distribution transformations, so any comment would be appreciated.

What sort of integration are you after? 2. Have you tried just using normalized counts?

The integration is supposed to give a better picture of transcriptomics as a whole, which would finally lead to specify the final destiny of a transcript based on its post-transcriptional regulation (whether it will translate or not). Till now we tried using TPM values to account for the different sequencing depths, as we have not used it yet to compare across conditions, but its behaviour across time. (Sorry, I erroneously flagged your comment while trying to upvote it :P)


text:
Differential expression analysis for a subset of transcripts
I want to check the differential expression of a specific class of transcripts (say, long non-coding RNAs) using DESeq2. Now, I know that the normalization step takes into account the total number of reads per sample. If I only input a subset of all the transcripts, the sequencing depth will change.
So my question is: can I simply input the subset of interest? Or should I input all transcripts, and subset the results later?

It's safer to import everything.  You'll want all the data for normalization and dispersion estimates.


text:
AutoDock 4 -soft docking
How could I turn on/off the smooth potential in AutoDock4? I need to do two docking. The first one is soft docking with smooth potential and the second is re-docking without the smooth potential.

The smooth parameter is an optional keyword and command that changes the smoothing potential. If not provided, the default is 0.5 angstroms by default. Try passing 0 with this keyword in the AutoGrid Parameter File to turn off smoothing potential.
See the AutoDock 4.2.5 documentation for more details.


text:
In silico PCR: any models that can quantitatively predict the initial content from the end amplified result?
Let's say I've run the PCR and I now have a certain amount of the amplicon product. Is there any model that infers how much of it was there before the PCR started? I suspect it should be a probabilistic model (a distribution) parameters (shape) of which will strongly depend on the details of the PCR protocol and the characteristics of the amplicon...
If not, then may be there are forward models that predict the amount of the amplified result from the initial content?

You need to know the primer efficiency for this to work, do you know that?

@DevonRyan Assume, yes. You can assume I know everything. At this point I want to simply familiarize myself with any such model if there are any.

Have you considered just using a qPCR calculator such as this http://nebiocalculator.neb.com./#!/qPCRGen ? The model of qPCR is pretty simple really, so there are a few of these premade.

This is qPCR based on CT number, it is common place to identify patient viraemia for example

The basic equation for PCR kinetics (Equation 1) states that the amount of amplicon after c cycles (Nc) is the starting concentration of the amplicon (N0) times the amplification efficiency (E) to the power c. The PCR efficiency in this equation is a number between 1 and 2 (2 indicates 100% efficiency).  NCBI
In this  article  on research gate they said  amplifx  was the software you were looking for to, calculate "(i.e. expected efficiency in PCR reaction)". The amplifx even mentions "Testing primers in in silico PCR".


text:
Searching for proteins by primary and secondary structure
I need to find peptides/proteins with a certain sequence of amino acids in a certain type of secondary structure. For example, "ALA PHE GLY" in alpha-helix. Is there any tool that does that?

Your request is rather specific, so there is no tool...
It is a straightforwards request though. I'd impose the requirement that the structures be solved and no predicted secondary structure. You can do it two ways if using Python.
PDB
After downloading all PDB structures, in python with PyMOL (installable with conda, tricky without) you'd do something like
import pymol2, os
from multiprocessing import Pool

def find_motif(file):
    motif = 'AFG'
    try:
        with pymol2.PyMOL() as pymol
            pymol.cmd.load(file)
            seq = pymol.cmd.get_fastastr()
            if motif not in seq:
                return None
            d = {x: []}
            # str.find will cause issues if there is a gap also structure may not start at one.
            pymol.cmd.iterate('resn ALA+PHE+GLY', 'x.append((resi, resn, ss))', space=d)
            sequence = -1
            i = 0
            for residue in d['x']:
                if residue['resn'] = 'ALA':
                    sequence = 0 # lazy enum
                    i = int(residue['resi'])
                elif sequence != -1 :
                    if int(residue['resi']) == i+1 and residue['resn'] = 'PHE':
                        sequence = 1
                    elif int(residue['resi']) == i+2 and residue['resn'] = 'GLY':
                        return (file, i)
                    else:
                        sequence = -1
    finally:
        return (file, None)

if __name__ == '__main__':
    p = Pool(12)
    structure_files = [f for f in os.listdir('pdb') if os.path.splitext(f)[1] == '.pdb']
    verdicts = p.map(finder, structure_files)
    print([(f,v) for f,v in verdicts if v is not None])

The end result will be a list of PDB with where in the sequence (in that PDB code).
(I wrote this in the reply box, so there may be a typo or two).
PyMOL has a selection algebra called pepseq, which matches missing atoms, but I was a whimp and played it safe by being more verbose. Also, a for loop can be used instead of an iterate. But I think it is slower.
atoms = pymol.cmd.get_model(f"pepseq {target} and name CA")
for atom in atoms.atom:
      if atom.ss == ' ...

Uniprot
In Uniprot, if a protein is crystallised, there will be the feature entries corresponding to the secondary structure. It is sheet or helix. No turns, pi-helices or 3_10 helices. There is no Biopython method for parsing Uniprot XML. I wrote myself a parser and it is not too horrid thing to do —although I monkeypatched elementtree. Alternatively, xmlschema.XMLSchema.to_dict can be used to convert the XML to a dictionary. The advantage of Uniprot is that the residues ids will not require correcting. The disadvantage is parsing XML.
My parser is specific for my application and it is poorly documented (so would not recommend it) but but... I can share it and human data. I'd be used like:
 import os
 from michelanglo_protein import ProteinCore
 ProteinCore.settings.startup(data_folder='protein-data')
 p = ProteinCore().gload(file=os.path.join('gpickle','Q86V25.gpz'))
 #print(p.asdict())
 if 'helix' in p.features:
        n = p.sequence.find( ...
        for ... p.features['helix'] #('helix', 'turn', 'strand')
        if x < n < y ...

Thank you. Parsing it on one's own doesn't look complicated. Still, I somehow balk at the idea that my needs are so specific... You are positive there's no such tool?

Proving a negative is hard, so I can only say I have never read of one.But I had a similar problem some time back and did it manually. Furthermore, logistically a server to do this would need its data pre-generated but then the combinatorial issues will be huge so would not work...

there is msdmotif https://www.ebi.ac.uk/pdbe/docs/Tutorials/msdmotif/introduction.html but I suspect the query might be a bit slow as there would be a lot of hits

I just realised that in my comment I said something incredibly daft. A server is totally possible, they just need to store sequence+SS type sequence as opposed to all combinations. Apologies.


text:
Normal score transformation
This question is related to a normal score transformation which is performed in the paper by
RUAN, Quansong, et al. Local similarity analysis reveals unique associations among marine bacterioplankton species and environmental factors. Bioinformatics, 2006, 22.20: 2532-2538.
They claim that given a variable X with observations $x= (x_1,\ldots,x_n)$ that are not normally distributed, normal score transformation is performed by first creating a rank vector observations $x$, $R^x = (R_1^x,\ldots, R_n^x)$ and then using and then transforming to $x'=(x_1',\ldots,x_n)$, by $x_i\ = \phi_{-1} (\frac{R_i^x}{n+1} )$ where $\phi_{-1}$ is the inverse of cummulative normal distribution.
It's intuitively clear (but maybe a bit harder to prove?) that resulting variable $X'$ has observations normally distributed, but I have a question about time series dependence. If we had that observations $x$ were correlated, say, over 4 adjacent observations, , i.e. $x_{n+1} = f(x_{n},x_{n-1},x_{n-2})$, how does this change after the normal transformation? Can we still find a function g, so that $x'_{n+1} = g(x'_{n},x'_{n-1},x'_{n-2})$ would hold? Since we're transforming based on the rank of data, I expect this to be not possible?

I'm voting to close this question as off-topic because is it general stats and beliongs in a different forum

This is biostatistics, the fact the was published in Bioinforamtics is an indication. The problem with the question is that none of the variables are defined so it is impossible to answer.

As you suspected, correlations will be heavily changed by such a transformation. The reason for this is also exactly as you described, namely that using ranks rather than values for the transformation eliminates all of the information needed to actually preserve these relationships. As an example, suppose we have a single observation with rank 1 across all observations that none the less shows a decrease in absolute value across them (i.e., there's a correlation with a negative slope over time). Such a transformation will result in a correlation of 0.
One might be able to find a different transformation that preserves this relationship, or might be able to simply pool all of the time points together before computing the ranks (this would then at least partially preserve some correlations, though depending on what's needed to properly normalize the data it may cause other problems), but in general there's no guarantee that this is possible.

It appears the original observations were not a timeseries and the OP is specifically asking about application to timeseries data.
In my experience these are two separate areas of analysis because a temporal skew may result in aggregate behaviour which violates the original transformation. For example hot weather versus cold weather.
What the OP appears to be doing is normalising each point in the time series, albeit I don't know what X represents. Personally I have doubts because the aggregate correlation is across all time points and one set of aggregrates, e.g. Ecoregion 1 is compared with another e.g. Ecoregion 2. This for example would be done via linear modelling and assessement of fit performed through many approaches, e.g. the residual is just one approach.
If data has large temporal skews your approach is unlikely to work, ultimately you don't know the variation in the time series which could far more heterogeneous than the agrregrate data. Seasonality is likely a major factor in your data and I would consider this in the first instance. There are lots of ways to transform a time series, standardisation is a classic approach. Again my personal experience is these are two separate areas of statistics and directly applying one to the other is something that should be assessed with care.

So just to try and explain in this scenario terms, when its cold if the variation in bacterioplankton is much greater than when its hot (smaller population sizes increases variance), i.e. at different times of year then your approach can't work because normalisation is about a fixed (constant) relationship between mean and standard deviation.


text:
Is New Tuxedo approach really better than the earlier ones?
I am planning to perform multiple routine NGS analyses for a hospital lab on patients' samples on an everyday basis.
Before I used a new tuxedo approach with ballgown etc and was rather satisfied with it.
I wonder if for the above-mentioned job it's ok to use a ballgown or, perhaps, that is not the best choice?

I'm fairly used to analyse NGS illumina data, but what is the New Tuxedo approach you mention? How do you measure best choice?

Thanks I am talking about this approach https://www.nature.com/articles/nprot.2016.095