text:
BWA doesn't recognize HiC options
I'm trying to run bwa on HiC data. This is the command I'm giving:
bwa mem -5SP -T0 -t4 -R '@RG\tID:\tSM:\tLB:\tPL:ILLUMINA\tPU:none' REFERENCE.fa PAIR1.fq PAIR2.fq

and here's the error message:
mem: illegal option -- 5

Supposedly, the 5SP option is for HiC data (cite), but my BWA isn't recognizing it. Anybody have any ideas? My initial thought is to update BWA to the newest version, but I don't want to mess with our replicatibility since I am new to the lab and don't want to be updating software when I don't need to. The server is a Mac Pro and bwa version is 0.7.12-r1039. Thanks.

You can use conda to install and test the desired version of bwa to test. This will not interfere with whatever version installed on the server.

This option has been added in 0.7.16a, see https://github.com/lh3/bwa/releases, so consider updating.

Thanks for providing exactly the answer needed!

Your link with those command line options is 1 year old.  Looks like that version of bwa is 6 years old.  You probably should update.


text:
Cobratoolbox unable to identify gurobi solver when passing initCobraToolbox
I am running into this issue lately when I am trying to install the Cobratoolbox on matlab along with having installed the prerequisite gurobi solver and adding it to the MATLAB path. I have run the gurobi_setup.m file as well. However I am unable to uninstall the cobratoolbox, as everytime i run initCobraToolbox.m it fails as it not able to recognize the solver. Is there anything I am missing in the installation? I have attached a screenshot for reference

In my opinion it is better to redirect the question to cobra-toolbox google group https://groups.google.com/g/cobra-toolbox

Not that I would know this package, but Installation should be clear in the README file and if you run into troubles while following it, it's good enough to open an issue on the platform where they develop the package. Not that I would think this question does not belong here, but I think you have much better chances there.

source: Stack Exchange

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

sort: votes


text:
How to get % similarity between strains and mutation files
I’m very new to python, and having some difficulty getting hang of some more complicated things
I have multiple files which look like so:

hCoV-19/Singapore/4/2020|EPI_ISL_410535|2020-02-03

hCoV-19/USA/WA13-UW9/2020|EPI_ISL_413601|2020-03-02

hCoV-19/USA/WA-UW142/2020|EPI_ISL_416680|2020-03-11

Please be aware that the lines above are meant to be one file
I want to extract the EPI_ISL_000000 for an easy comparison among files.
Could someone please advise on:

A programme to extract this data into new files (There’s many lines in each file - 1000+)

A programme to then give a % comparison between two or more files - comparing all lines in one file against all lines in a second+ file

Anyway I get it, the first part is reg-ex, either via Python's re module, sed or a perl onliner. The second part, is easier to do in MEGA X.

I have now figured this out, thank you!

If you have figured it out please answer your own question and accept it. thank you.

I didn't study this subject... however, the percentage between two organisms may depend heavily on local protein aligniments and the same other active genetic material (because not active rna or not protein make rna has to have weight 0)

compile proteins of virus1 compile proteins of virus2 search % equality of proteins of the 2 viruses... or align the genome of 2 viruses and find uguagliance i don't know what number is more significative

   left_lineagelist = [x.split('_')[-1].split('|')[0]
                          for x in left_lineagelist]
        right_lineagelist = set([x.split('_')[-1].split('|')[0]
                          for x in right_lineagelist])

Allows for extraction of 6 digit EPI, provided the file has had sequences removed prior; as such:
for line in lines:
    if line[0] == '>':
        print(line[1:])

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


text:
dbnSNP frequency anomalies
Sometimes dbSNP reports very different allele frequencies for different large-scale genome projects e.g. between 1000 Genomes and GnomAD
rs11822440 1000Genomes A=0.4629    C=0.5371 GnomAD A=0.99997   T=0.00003
rs28531832 1000Genomes T=0.4379    C=0.5621 GnomAD T=0.99993   A=0.00007
rs78931618 1000Genomes A=0.7334    C=0.2666 GnomAD A=0.9999    T=0.0001
rs75716913 1000Genomes T=0.9006    C=0.0994 GnomAD T=0.99149   A=0.00851
Is this an issue with one or other of these projects data processing? Is something else going on here biologically or technically that might cause this?

do you know what populations were used to calculate the gnomad frequencies?

These are SNP frequencies - its not the same thing as whole genome GC content. What we are looking at here is selection pressure at given SNPs

@user438383 the pattern is consistently different between projects in global and sub-populations e.g. African, European

@Michael these are the same SNPs with different frequencies for different projects - but both on "normal human" so I would be surprised if there is that much selection via recruitment.

The other thing - have you double checked you are comparing the same genome build between the two (i.e. gnomad is gr38 and 1kg is gr37)?

allopatry .. it happens

I'll stick my neck out and simply state that allopatry alone would account for the pattern you are observing. Thus populations can exihibit variation in allelic profiles, particularly where Mendellian inheritance is operating (i.e. higher eukaryotes), singly due to geographic isolation.
This represents the Meyers view of population genetics and would account for variation within neutral theory, e.g. silent mutations. Obviously I'm guessing the populations were selected by geographic foci in generating the allelic frequencies.


text:
Trimmomatic trims most of the reverse strand
I have a paired-end data. It is from small part of human genome not a WGS or WES. I use BWA for alignment and do variant calling. Since I had some false positives I wanted to do trimming.
I used trimmomatic with command:
trimmomatic PE -phred33 -threads 2 read1.fastq.gz read2.fastq.gz read1_trimmed.fastq read2_trimmed.fastq LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15`

I do not know why but I have this output:
Input Read Pairs: 1834513
Both Surviving: 1822034 (99.32%)
Forward Only Surviving: 5740 (0.31%)
Reverse Only Surviving: 417 (0.02%)
Dropped: 6322 (0.34%)
TrimmomaticPE: Completed successfully

I do not have any adapters since the are removed in bcl2fastq step. In FASTQ quality control I have seen that quality of reverse and forward reads are similar in quality but trimmomatic trimmed many reads in reverse file. I cannot do alignment with these files, BWA gives a warning that says:
[bseq_read] the 2nd file has fewer sequences.
[process] read 155012 sequences (16693324 bp)...
[bseq_read] the 2nd file has fewer sequences.

Is there anyone who faced a problem like this before?
Thank you in advance

Maybe you need to remove the unpaired reads:
trimmomatic PE -phred33 -threads 2 \
    read1.fastq.gz read2.fastq.gz \
    read1_trimmed.fastq.gz read1_unique.fastq.gz \
    read2_trimmed.fastq.gz read2_unique.fastq.gz \
    LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15

Then you can try with bwa.


text:
WGBS vs ONT, which one should I trust?
First of all let me apologize if this is not the place for this question. While it can be quite broad, I will try to make it as specific as I can.
Context
I am testing the correlation between Whole Genome Bisulfite Sequencing and Oxford Nanopore methylation calls for one individual (mammal, blood tissue) with replicates at all CpG sites. In the (very) best cases I find a close to 0.9 correlation. Yet, for the part that does not correlate I would like to know which technology is closer to the "biological truth".
My strategy
My plan would be to first find CpG sites/short CpG regions showing differential methylation states between WGBS and ONT. Then I would validate the methylation states of these sites using a third party technology such as pyromark.
What I am asking to the SE community

First, has anyone done this before ? I would love to hear about your experience

Do you see some critical aspects I should be watching out for ?

Last if you have comments about the strategy I propose here, you are welcome

I'd probably add EM-Seq (https://www.neb.com/products/e7120-nebnext-enzymatic-methyl-seq-kit) to the comparison, since it appears to become a valid replacement for WGBS based on enzymes rather than bisulfilte conversion.

@SebastianMüller Thank youfor the suggestion, I didn't know about EM-Seq so I will look at it.

Methylation can be cell-specific, which makes it difficult to evaluate accuracy on a bulk-cell level (even within the same tissue). How can you tell that the differences you're seeing are due to platform differences, or due to biological variation?
I find that adding more haystacks doesn't help much in working out the truth of a dataset. If you want to investigate biological truth, it would be better to create a biological system with a known methylation (or demethylation) pattern and test that.

thank you for your interest. Actually these datasets come from the same tissue sample. But you are right, working on a bulk-cell level makes it difficult to distinguish tech bias from biological variation. I think it can be difficult to produce control data that fits the biological variation you find in a bulk of wild cell, This is why I thought using a third party technology.

I would assess directionality and accuracy of prediction by 1) WGBS predicting ONT and then 2) ONT predicting WGBS.
Firstly, I would use deep learning (or machine learning) and train WGBS against ONT, parameterise and then test. Then conversely train ONT against WGBS, parameterise and then test. The approach with highest accuracy of prediction (using the 'accuracy' index) would be the approach assumed to be "closer to the biological truth". If both calculations produced comparable accuracies I would conclude what @gringer has stated that natural variation / natural heterogeneity is the predominant signal in the signal. This conclusion depends on the ability of deep learning to map the biological process.
It is an approach highly applicable to a deep learning estimate regardless and may not replace a truely controlled test, this will depend on what controls you have run, but could provide a valuable and easily obtainable insight in its absence.
If you did run a control sample, then you are truely in, because this would provide the training set and the WBGS and ONT provide the test data.

You asked me how it works
It doesn't replace good controls, but it is 'trendy workaround'.
Two steps

In deep learning (artificial neural networks, ANN, in this case you would use RNN such as LSTM) or machine learning (ML) e.g. random forest - random forests are MUCH easier to do, you supply a series of known values, thus CpG sites and their WGBS scores and use ANN (RNN approach) or ML (e.g. random forest). Both methods will look for pattens between the observed DNA data and the CpG occurence (training set). Usually this is 60% of your data set is used in training. A separate 20% of the data set is often used for parameterisation and a final 20% of the data set is used for testing. Thus how well can a ANN or ML predict a CpG site based on WGBS data. Likewise you would do the same for ONT. There is a caveat here which I discuss below.

The idea is if you give it a bit of DNA it will assign whether it is CpG (see caveat below).

However, instead of testing WGBS to WGBS you also test it WGBS to ONT and ONT to WGBS and the method with the highest 'accuracy' (positive hits, against true hits) is 'better method'. This is because you already know which is ONT positive and ask a WGBS trained algorithm whether it thinks its positive. The relationship between observed CpG and 'true (experimental) CpG is your accuracy score.

The caveat is that you need to supply the training set with sites that are known to be negative under WGBS or ONT. You need to think the best strategy for doing this. If some of those negatives were positve in the other method that would affect the outcome alot.
There is an issue about vectorisation this is often done using k-mers, and is a headache. Basically you have to make a bit of DNA into numbers and that in my opinion is the difficult part. In your case I ain't so sure its so difficult. If you can get that bit correct and biologically meaningful you are in business. My stuff is phylogeny based and I understand the relationship between biology (mutation) and numbers. In your case that needs thought, however someone will have thought about it and solved it.
As a personal statement I would do this via random forests, because ANN is difficult.
One such vectorisation method is here by Li et al (2017), however as they have used it, the method has flaws due to equal weighting of mutation frequencies. However, these flaws do not apply to your work - it would actually be a good method, they singly apply to tree building.

The final thing is that the more experimental data you can feed it the better. You don't have to complicated about it, just create a new column in the input section, e.g. along side teh vectorised DNA and the algorithm will figure out whether it helps. You can't feed it enough data, because it is doesn't help it wouldn't use it.

Thank you Michael, I am really interested by your approach but I think I don't understand what you mean by predicting WGBS (though I know the rôle of deep learning for ONT, I trained some models). I also have control samples (but the fully methylated one may not be trustable) .

@PaulEndymion I've modified the post above


text:
How to find Genotype, from list of mutations:lineage list?
I'm very new to python and I'm getting along nicely (I'd like to believe); however, I must be missing something here. I'm looking to read each and every file and compare each line against one another.
The aim of this is to see if any mutations are shared in any of the strains - to determine the genotype.
The inner-workings of the files are single spaced headers, as those you would find on GISAID.
os.chdir("C:)
with open('A D614G.txt', 'r') as file2:
with open('A E780Q.txt', 'r') as file3:
with open('A G476S.txt', 'r') as file4:
with open('A L18F.txt', 'r') as file5:
with open('A N439K.txt', 'r') as file6:
with open('A S477N.txt', 'r') as file7:
with open('A T478I.txt', 'r') as file8:
with open('A V483A.txt', 'r') as file9:

files = [file1, file2, file3,file4, file5, file6, file7, file8, file9]

it = itertools.permutations((files), len(files))

for x in it:
    print(x)

os.chdir("C:)
with open('tester.txt', 'wt') as file_out:
    for x in it:
        file_out.write()```

1 I recieve the following error on a very long repeated stretch of text

```<_io.TextIOWrapper name='A G476S.txt' mode='r' encoding='cp1252'>```

I seem to have worked it out - I shall post the working code for anyone interested. If anyone would like to streamline or test for efficiency please feel free. Very proud of myself

If you've answered your own question, you should add it as a an answer and then accept it so that other people know it works.

I do agree with this suggestion of posting your solution in the 'answers' section @TheoJones

I’ve been looking for the code, it’s deep in my files somewhere. It’s also the most horrible looking thing, so be prepared

import re
import os
import pandas as pd
import logging
import sys
import csv
import itertools

mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

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

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

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

combination_labels = []
combination_counts = []

for lineage in lineages:

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

    for combination in combinations:

            combination_df = df[list(df)]

            for mutation in combination:
                combination_df = combination_df[combination_df[mutation] == 1]

        #print(combination_df.shape[0])

            combination_labels.append('_'.join(combination))
            combination_counts.append(combination_df.shape[0])

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

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


text:
Performing PCA for the samples and for the genes
I have 10 samples from a RNAseq experiment (5 control, 5 disease), I have performed a cluster analysis for the samples and for the genes (4000 genes aprox) to see how they cluster (to see which samples are similar and which genes have similar expressions).
I was wondering if in terms of statistics it makes any sense to perform a PCA of the samples instead of the genes?
thank you in advance!

I think doing PCA on RNAseq samples is something people do - you can probably just use the prcomp function in R. However, I'm not sure how much information you can get in a PCA with just 10 samples.

People usually do PCA on RNAseq samples, but they use genes as variables. But I want to do it using samples as variables, and I am not sure if that makes sense in terms of statistics.

Why do you specifically want to do a PCA using samples? Can you do it on genes and plot the loading of samples? Since you have thousands of genes and only 10 samples, it doesn't seem you can do it on samples.

First of all, I am not a statistician but have been performing RNA-seq analysis for a while, here is my take:

I think performing PCA on the samples makes sense in terms of mathematics. What I mean is you can compute eigenvectors and subsequently principal components either on your genes or on your samples.

I don't think performing PCA on the samples makes sense in terms of statistics in RNA-seq context (and probably in any other context): Your genes (features) describe your samples and not the other way around. So the computed PCs should describe/explain your samples and not the other way around. You would not perform PCA on your samples (individual cells in scRNA-seq) even in the case of scRNA-seq, where number of samples sometimes exceeds the number of genes, there PCA is performed on genes and the resulting few PCs are used for downstream analysis, usually clustering and visualization, both used to infer about samples (cells).

If you use, say, R's prcomp, you can use my_pca_object$x to get PC coordinates for the samples, and my_pca_object$rotation to get how much each gene contributes to each PCA.  (You might need to transpose your transformed count file to get it to work right)

It is quite possible to make a PCA plot with the samples. Follow this tutorial : https://hbctraining.github.io/DGE_workshop/lessons/03_DGE_QC_analysis.html


text:
.cif .pdb conversion with python
can you please give me an advice how to convert .cif files into .pdb preferably using python?
Thanks in advance!
Best,
Balint Biro

Its certainly doable and PyMol will also do this. A number of site members hold this expertise. What I think is that the format you are downloading in is cif rather than what you want which is pdb. If that is the case please modify your question

Dear Michael,yes, that is the case but I cant see how should I modify my question...

Using the [edit] button :-)

As @Michael said PyMOL can do this. PyMOL is not only an app, but also a python package —installed via conda not via the regular download (you can have both).
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load(infile,'myprotein')
     pymol.cmd.save(infile.replace('.cif', '.pdb'), selection='myprotein')

This runs the parallelisable instance, so you can have that context manager block running on how many threads or processes you fancy.

Dear Matteo, thanks for that. Unfortunately it gives me back the error that 'name 'infile' is not defined'.

@BalintBiro You need to change that part so that it matches the name of the file you are using.

Hi @BalintBiro please upvote or mark 'accepted'. This is the answer. You need to define the location of the infile before the with statement starts, thats why you have the error. You need a certain degree of coding proficiency.


text:
tissue cross talk using gene expression data
I am studying about tissue cross-talk using gene expression data. My first goal is to identify the genes responsible for communication signaling(cross talk).  Let say the tissues are fat and muscle.
One thing that came to my mind was doing correlation analysis and visualizing using a heatmap.
Here I would like to get your help with the following points:
How to extract the one where both tissues showing a possible correlation from the four possible quadrants of the correlation matrix?
How best I can visualize that quadrant using heatmap or any other plotting methods.
Here I simulated the data as follows. Rows are genes and columns as samples.
set.seed(1)
A <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(A) <- paste0("G_", 1:100)
colnames(A) <- paste0("S_", 1:5)

set.seed(42)
B <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(B) <- paste0("G_", 1:100)
colnames(B) <- paste0("S_", 1:5)

Tank you!

Is there any evidence that gene expression can inform about crosstalk? I mean we are talking about some kind of ligand-receptor interaction here, or receptor-receptor interaction. That is the realm of proteins and protein modifications, signaling cascades etc. I am not sure that looking at the transcriptome will give you any information here.

@ATpoint; I know it is tricky. I plan to do proteomics, but this is to get an idea for proteomics experiments.


text:
Power law distribution alpha values
I did a Pearson co-expression analysis for generating networks for my tissue-specific (chondrocytes) RNA-seq data. I used R package poweRlaw to check for the power-law distribution. We got the following values:
 Xmin = 14
 goodness of fit = 0.035
 α = 4.76

My question is what does alpha (α) of more than 4 indicates?
Is it is okay to have such values?

Thank you

It means the distribution is heavily skewed and has a 'long-tail', thus the variance is much greater than the mean. In other words most values are around zero, with a small number of values having extreme values.
[
You need to check the goodness of fit because that is saying the deviation from a power-laws distribution is just significant. So it ain't quite power-laws, but the general gist of a 'long-tailed distribution' is correct.
Its a good result and a nice analysis, but you are not quite there at present (needs to be >0.05). Check the negative binomial distribution.

I've seen you data, it looks very good and the power series is not a bad explanation of the data. There are two possibilities for the goodness of fit being 0.035

The goodness of fit is low because of 'outliers' at the beginning of the distribution
An alternative distribution, such as a negative binomial distribution fits better
Check it is 0.035 not 0.35 because that changes everything.

The best way to do this is to plot the power series expected distribution onto the graph and look where the fit 'fails'.
If it is not fitting at the beginning of the distribution thats is about outliers and easily solvable. If it is not really capturing the shape then its probably a slightly different distribution.
If its point 1 then just remove the very first point in the graph and try again and check the p-value. Or take the mean of the clump of three points next to the first point and try again (or just remove those three point). Don't do any more than that and the technical word is "outlier" analysis. In your write-up you must explain the issue around any 'outliers' and that is perfectly acceptable. If it is a power-series distrupted by an outlier(s) there is formal analysis to demonstrate an 'outlier' using outlier analysis if any one ever asks you.
If its point 2. then use e.g. negative bionomial or the gamma distribution (alpha <1) and check the goodness of fit.
You want a goodness of fit > 0.05 and preferably >0.1 to be sure its the right distribution. There are formal statistics to do this, but that is the main one.
My guess is there are outliers at the beginning because your sample size is not large enough to identify the variation correctly, hence the fit is low. In 'long-tails' its easy to have a 'low sample size' because of the way they work. Thus in this case mapping a distribution is a perfect way to compensate for the sampling.

Thank you for your reply. I have now attached my result graph and I would really appreciate it if you could have a look at it

Hi @Priya I've modificed the answer above

Thank you very much for explaining it so nicely.

@Michael, it is not clear to me whether alpha here represents skewness or something else. Similarly, goodness of fit interpretation does not seem to help.

The power laws distribution is skewed, but alpha will be 'more skewed' as the value increases. The value you have will be optimal for your the data set you plotted. I'm not sure what you mean by goodness of fit "doesn't seem to help".


text:
Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper
I am trying to replicate the results of this paper roughly guided by this pipeline.
Basically, we are trying to detect a differential expression between mesenchymal and epithelial cells under cisplatin treatment vs control.
However, I am facing a problem when trying to calculate the p-value for the Mann-Whitney U-test. This is the part of the paper I am referring to:

Robust Multichip Average (RMA) normalization was performed at the
transcript level on the results from the Affymetrix Human Gene ST1.0
arrays using Affymetrix Power Tool 1.14.4.1 for all 46 sham- or
cisplatin-treated ovarian cancer cell lines. The normalized data were
subsequently standardized using ComBat71 to remove the batch effect.
In this experiment, the cisplatin treatment assay was performed in
triplicate on 20 cell lines, while single assays (without replicate)
were performed on the remaining 26 cell lines. Taking advantage of the
triplicate data, potentially fragile probes with strong variations (an
s.d. of >0.2) within the triplicates were removed, decreasing the
probe number from 33 297 to 21 329. To perform a fair comparison, the
triplicate data were then log-averaged into one value so that one
result for each cell line could be used in the following analyses. To
detect differential responses to cisplatin between epithelial- and
mesenchymal-like cell lines, the transcriptomic responses to cisplatin
were computed by subtracting the gene expression value of control
(cisplatin untreated) cells from that of cisplatin-treated cells.
Mann–Whitney U-test (P<0.01 as a cut-off value) was subsequently used
to detect the differential transcriptomic responses between the
expression changes by cisplatin treatment in epithelial-like cell
lines with those in mesenchymal-like cell lines (Supplementary Table
3).

My p-values for a given gene are astronomically different to that of the paper...
Here is my code:
expr_f is the expression values matrix from the Expression set (log-transformed).
SDRF_f is the pData from the same expression set.
I am feeding into the wilcox.test function a vector of expression from epithelial cells (epi_values) and mesencymal cells (mes_values) for the same probe (8143663).
 genelist=rownames(expr_f)
epi_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Epithelial-like"))
mes_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Mesenchymal-like"))

epi_values = expr_f["8143663", epi_samples]
mes_values = expr_f["8143663", mes_samples]

wilcox.test(mes_values, epi_values, exact=TRUE)

In the paper, their p-value for this probe is 0.000372259258165627 while mine is 0.8785. The difference makes no sense at all to me.
I don't know if I am not quite getting the concept of the u-test, if I am using the function wrong, or if it's something else. Any input is appreciated

Cross-posted on biostars: https://www.biostars.org/p/467695/

Wow that is some difference particularly for the Mann-Whitney. You are using Wilcox of course. Regardless there is one possibility which I will mention below if I get chance.

I think you may need to check your pipeline in this regard. I ran the analysis and replicated as closely as I could their M&M (I didn't filter the "fragile probes" since you anyway care only about a specific probe). In any case, I get something relatively close to their result (P=0.004059), which is still one order of magnitude off, but without more detailed M&M, this is probably as close as I can get to their analysis. Here's what I did (if you prefer to read the notebook html, that is here)
library(Biobase)
library(GEOquery)
library(sva) # ComBat
library(oligo)
library(googlesheets4)
library(stringr)
library(dplyr)

#' Get the GEO dataset mostly for the phenoData annotation
gset <- getGEO("GSE47856", GSEMatrix =TRUE, getGPL=FALSE)
adf <- Biobase::phenoData(gset[[1]])
sampleNames(adf)

#' Read RAW cel files from the GEO data supplementary rar archive
#' https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE47856&format=file
#' In my case extracted to ~/Downloads/cel
#' Ensure that we read them in the same order as the phenoData object
celfiles <- sapply(sampleNames(adf), function(f) {
  pat <- paste0(f, ".*.CEL.gz$")
  dir("~/Downloads/cel", pattern = pat, full.names = TRUE)
})
raw_data <- oligo::read.celfiles(celfiles,
                                 phenoData = adf,
                                 experimentData = experimentData(gset[[1]]))

#' Get the raw expression and RMA on transcript clusters
rma_expr <- oligo::rma(raw_data, target = "core")
rma_expr_assay <- assayData(rma_expr)$exprs

#' Use ComBat for batch effect removal
batch_data <- phenoData(rma_expr)[["batch:ch1"]]
cb_expr <- sva::ComBat(rma_expr_assay, batch_data)

#' Get supplementary table 1 from the paper (which I copied on this Google Sheet)
supp_data <- read_sheet("https://docs.google.com/spreadsheets/d/142__1O35zvlttejvz62XQX33uLRJhqIDvb604smVSRo")

#' Match the 'cell line' phenotype from the phenoData to the sample info of
#' the supplementary table and extract the cell classification. Also remove
#' all white space and convert everything to lower case for easier matching.
#' Some cell lines also have a/b attached, so I'm dropping that to match them
#' up. Mayb not be correct to do that here and treat all a/b the same (?)
cell_line <- str_replace(
  str_replace(
  str_remove(
    str_to_lower(adf[["cell line:ch1"]]),
    " +"),
  "^tay$", "taya"),
  "dov13[ab]", "dov13"
)

sample_info <- str_remove(str_to_lower(supp_data$Sample), " +")
classification <- supp_data$Classification[match(cell_line, sample_info)]

#' Get the treament from the phenoData object
treatment <- adf[["treatment:ch1"]]

#' We'll need the geometric mean later (R has no built-in)
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

#' Set up a tibble (==data frame) with all that info for your probe
tibble(Expression = cb_expr["8143663", ],
       Line = cell_line,
       Treatment = treatment,
       Classification = classification) %>%
  # Summarise replicated data with the geometric mean
  group_by(Line, Treatment, Classification) %>%
  summarise(Emean = gm_mean(Expression)) %>%
  # Make sure data is in order
  arrange(Line, Treatment, Classification)%>%
  # Calculate the difference between the cell types
  group_by(Line, Classification) %>%
  summarise(EmeanDiff = diff(Emean)) -> xdata

#' Test the difference
xdata <- split(xdata$EmeanDiff, xdata$Classification)
wilcox.test(x = xdata$`Epithelial-like`, y = xdata$`Mesenchymal-like`,
            exact = TRUE)

#' Session info
sessionInfo()
```

Hi @CSP it is certainly an interesting discussion and some good R coding. Could you look through this very detailed reply and comment please?

You were right, I had a minor mistake in my pipeline! Thank you for this answer, very detailed and useful!

Great that it worked. If this is answered, please accept the answer

Agreed @CSP please accept the above answer

The most likely explanation here is paired versus unpaired.
Mann-Whitney is normally associated with unpaired data, but the authors may have performed a paired test. You performed an unpaired test and thus there is a massive discrepancy in results. The authors should have stated the effect of pairing - it that is what they did - but they are under no strict obligations to do so. You however have performed a Wilcoxan rank sign test, albeit the results shouldn't be so different.

In R, a two-sample wilcox.test (i.e. providing both x and y) is the same as a Mann-Whitney U test. From the docs:

Performs one- and two-sample Wilcoxon tests on vectors of data; the latter is also known as ‘Mann-Whitney’ test.

Also they cannot do a paired test, since they have uneven groups (Epithelial vs Mesenchymal)

Okay fair point @BastianSchiffthaler , I'm wrong. I was aware of the ability of Wilcox to approximate Mann-Whitney ... if the samples are naturally unpaired ... again simply I'm wrong. Mann-Whitney is naturally unpaired so I was never really sure. The only thing I'd point out is that generally any bug can result in p<<0.01 but when p>>0.1 it difficult for a bug to achieve this.

The samples are unpaired, it was a mistake in my prior pipeline, but thank you for you very useful input!


text:
PacBio long-reads impact in transcriptome de novo assembly?
We are strongly interested in assembly a good transcriptome of reference for a non-model organism and build a local database. We have sequenced the same individual with Illumina (150 millions of pair-end reads) and PacBios IsoSeq v3 (2 SMRT cells, one for shorter transcripts, shorter than 5kb and other for longer transcripts, up to 5kb).
To process long-reads, I have followed the PacBio IsoSeq pipeline proposed in their Github repo (https://github.com/PacificBiosciences/IsoSeq). The final result was removing 70% of the long-reads. Is that normal?
Using this data, I have assembled the transcriptome using only short reads and another combining long- and short-reads. In the end, I have not found any difference... Approx. The same N50, the same number of transcripts assembled, the rate of misassemblies... Does anyone know if PacBio data does not worth for transcriptome de novo assembly?

Could you clarify your question please ... " Does anyone know if PacBio data does not worth for transcriptome de novo assembly?" I think there is a spelling error which confuses the meaning.

Biostars cross-posting.

When are you saying removing, in which stage it got removed?

A few comments:
Never use N50 as a metric especially for transcriptomes. It has some semblance of relevance for genome assembly, but all that is void for a transcriptome with inherently dynamic lengths.
At the end of your IsoSeq pipeline, you should have (ideally) full-length transcripts. Have you considered that you don't miss much in your IsoSeq data? The biggest gain using IsoSeq is that it should have fewer non-existent isoforms assembled than a De-Bruin graph approach.
Let's say you have the following locus in your genome:
5'                                                                 3'
------------=======-----------=======----------=======-------------
Upstream      E1     Intron     E2    Intron     E3      Downstream

There are theoretically 7 isoforms:

E1
E2
E3
E1,E2
E2,E3
E1,E3
E1,E2,E3

A short read assembler has to look at evidence that connects exons to see which isoforms are represented. Let's say you have the following PE reads:
5'                                                                 3'
------------=======-----------=======----------=======-------------
               -->..............<--- Read 1
                                 -->...............<-- Read 2
                 --___________-->..<-__________--- Read 3

Read 1 and 2 both have a mate in either exon, so that's good evidence that they come from a single transcript. You can also have reads like Read 3, which spans two exon-exon junctions, again good evidence. There are some problems, however, if things get more complicated: if the cumulative length of your exons is greater than the fragment length of your short reads, or you have many exons in a gene, you will never get a single read pair that confirms the existence of some transcripts. The same goes for the opposite. If you find a read pair on a single exon, that is not evidence that this single exon is a transcript in and of itself.
Here's where short read assemblers take different approaches. Some make educated guesses about the existence of some transcripts, e.g. by utilizing expression counts in an EM or dynamic programming approach, some just enumerate possibilities.
And that's where IsoSeq comes in. A single, long read from a single transcript helps a lot to add evidence and make sure your transcripts actually represent the transcripts that exist in your organism. So, often, your IsoSeq should be fewer transcripts than your short read assembly.
In essence, what you should do, is take the IsoSeq and curate those as a set of high-confidence transcripts, i.e. those that you have observed from a single read source. That's not to say that some transcripts that exist only in the short read assembly aren't real, just that you have less definitive evidence for them.
If you are mostly interested in the gene level, that doesn't make a huge difference, but once you need to look into transcript level usage/expression etc it helps a lot to be able to eliminate low confidence transcripts.
Use tools better suited for transcriptome assembly evaluation:
Detonate, TransRate and (Busco) and compare your assemblies using those.

Hi! Thanks, I know that N50 is not the best proxy to evaluate my assembly. I've used BUSCO, transrate and also E90N50 as proxy to evaluate the differente assemblies and the result with and without long-reads is the same. may it be possible assembled all the isoforms using only short-reads?

That is entirely possible! I edited my answer to highlight the real advantage of isoseq data in more detail, as imo its greatest strength is that you have evidence for the existence of transcripts from a single-read source. It's a metric that is hard to benchmark so I'm not aware of any tools to assess that besides visual checks


text:
Apply trained PRS on another dataset
I am using PRSice to compute the PRS over a train set and want to use the coefficient used on the train set to apply it on another set which I will call the test set.
Once I compute the PRS I get a file called analysis.prsice
This is the few first lines

Pheno Set Threshold R2 P Coefficient Standard.Error Num_SNP

- Base 5e-08 0.00550157 0.0659631 -21.0513 11.4492 43

- Base 5.005e-05 0.00396637 0.113759 -61.8545 39.1108 184

- Base 0.00010005 0.00236122 0.220463 -68.9527 56.2742 286

- Base 0.00015005 0.00113035 0.394503 -59.4058 69.7677 382

- Base 0.00020005 0.00072253 0.495129 -55.2605 81.0064 471

- Base 0.00025005 0.000334179 0.642121 -41.5674 89.4432 538

- Base 0.00030005 7.5306e-05 0.825111 -21.391 96.8021 601

- Base 0.00035005 3.95773e-06 0.959527 5.25804 103.612 668

- Base 0.00040005 0.000175506 0.735018 37.723 111.455 743

I thought that I could apply the PRS on another dataset by multiplying the coefficient in the result by the allelic dosage of the corresponding SNPs. However, I am not sure how to know which are the corresponding SNPs of the coefficients.
Also, if there was a command that does that and it would be easier, could you please tell me how to do it?
Thanks

How does PRSice work, i.e. what ML is it using? What is the accurarcy of your predictions? Applying a trained data set to predict another data set should be straight forward. You don't need to specify individual coefficients

I am not sure what you mean? actually, I am not even sure how you measure PRS accuracy. I know in the C+T method that you add the allelic dosage weighted by the effect size in the GWAS

The coefficient is the third column from the right, you can tell because some are positive and other negative. When the value is negative the association is a negative correlation. When it is positive (last two) that is more like the result you are seeking. The further away from zero the stronger the correlation (positive or negative)
These are very large values for the coefficients. The SNP is the last column on the right. I would disgard the top two results because the probabilities are too borderline, e.g. if you used Bonferoni's correction.
There isn't enough information on this calculation nor the data set to make any further insights other than you have 7 SNPs showing PRS


text:
Difference between paired-end, mate-pair and long read
I writing here because I have some questions for you.
I wondered what the essential differences were between paired-end, mate-pair and long read?
for me mate-pair and classic paired-end are both paired-end reads, with the difference that :

For classical paired-end: the insert size on classic paired-end is
smaller (about 500bp)
while the insert size of mate-pair is much longer (several Kb)
which allows to join the contiguous between them especially is it?

In revenge for the long-reads, I imagine that they are simply reads that are synthesized with a large read size but that do not allow like the maite-pair to make junctions between contigs?

https://www.biostars.org/p/467845/

They are all very different in separate regards, but they all refer to different wet-lab and sequencing protocols/technologies.
First, PE (paired end) reads are typically short (50-300) reads, most often Illumina HiSeq, MiSeq or NovaSeq protocols. Both pairs originate from a single fragment which is sequenced from either end:
------> PE1
-------------------------------- Fragment
                         <------ PE2
       ^^^^^^^^^^^^^^^^^^        Insert

In contrast mate pairs arise from a fragment that is circularized before sequencing. Now imagine the fragment is circularized:
                     -----> MP1
....--------------------------------.......
          MP2 <-----
                    ^Circularized (joined) here

Long reads typically refer to PacBio or Oxford nanopore data, with very long reads lengths (up to several hundred kb).

thank you very much for your time, as you can see ATpoint gave a link where I also asked a question in respond to a post, sorry for the duplication...


text:
Please help me with my command!
I am looking for a command to count the total number of unique proteins in a file. For instance I wrote a command to BLAST XYZ proteins as query against the DLY proteins as a database to determine the number of genes or proteins that are shared between XYZ and DLY:
blastp -query ~/user/XYZ.faa -db DLY -evalue 1e-10 -outfmt '6 qseqid' | sort | uniq | wc -l

I used the -evalue option to provide me with alignments with an e-value that is 10^-10 or lower.
How can I represent my command so that that it can provide me with a proportion of the total proteins in each species (i.e. XYZ and DLY).
Thank you!

I notice you changed your title from "how to X" to "please help with X!" actually the former is better, however, your title is not informative at all hence the very low view counts (it is the weekend though). "Tallying the proportion of the total proteins in each species from a Blast search" would be a more informative title. I am human so I do things here because they are interesting —so you need to convince me an other volunteer to look at your post. Same applies with Twitter posts, papers etc. informativeness/catchiness gets you engagement.

Anyway, about your actual question. It is best done programmatically. Species are in NCBI in square brackets in the header. So grep with -o flag and a regex pattern say '\[.*?\] would get those brakets, as you don't want the square brackets, a lookbehind will probably do, e.g. grep -o "(?<=\[)[^\]]+" -. However, there are several problems with a species tally of a blast result...

If one wanted the count of paralogues, the sequenced species are not a reflection of diversity —if you search anything E. coli you get thousands of strains. In your case, orthologue count per species, has the problem that complete genomes only would work and even then you'd have to make sure of this... So sed or grep are a poor choice. Doing it manually for representative species is actually the easiest solution...


text:
Error in cor(exprs(gse), use = "c") : no complete element pairs
I'm trying to plot a heat map for a "GSE133399", I retrieved the GSE using GEOquery package and then assigned to an object and then tried to plot a heatmap, but I was not successful, I used the following code in RStudio:
install.packages("BiocManager")
BiocManager::install("GEOquery")
library(GEOquery)
## change my_id to be the dataset that you want.
my_id <- "GSE133399"
gse <- getGEO(my_id)
length(gse)
gse <- gse[[1]]
gse
pData(gse) ## print the sample information
fData(gse) ## print the gene annotation
exprs(gse) ## print the expression data

## exprs get the expression levels as a data frame and get the distribution
summary(exprs(gse))
exprs(gse) <- log2(exprs(gse))
boxplot(exprs(gse),outline=FALSE)
library(dplyr)
sampleInfo <- pData(gse)
sampleInfo
## source_name_ch1 and characteristics_ch1.1 seem to contain factors we might need for the analysis. Let's pick just those columns

sampleInfo <- select(sampleInfo, source_name_ch1,characteristics_ch1.1)

## Optionally, rename to more convenient column names
sampleInfo <- rename(sampleInfo,group = source_name_ch1, patient=characteristics_ch1.1)
sampleInfo
install.packages("pheatmap")
library(pheatmap)
## argument use="c" stops an error if there are any missing data points

corMatrix <- cor(exprs(gse),use="c")

I get this error:
Error in cor(exprs(gse), use = "c") : no complete element pairs

I also  tried this:
corMatrix <- cor(exprs(gse),use="c", use="complete.obs")

I get this error:
Error in cor(exprs(gse), use = "c", use = "complete.obs") :
  formal argument "use" matched by multiple actual arguments

I want to achieve this plot:

The GSE133399 dataset contains both ATAC-seq and RNA-seq data, none of which can be queried with GEOquery which works only for microarray data.
exprs(gse) ## print the expression data returns an empty object as you might have noticed, therefore the errors in the lines that follow this command.
On the bottom of https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133399 at the Supplementary file section you can download FPKM-transformed count data and use these. I would be careful though as FPKM is not ideal for inter-sample comparison. I personally always (if feasable) download the raw data and then process myself, but for a guess on what is going on the FPKMs probably will do.

How do you process the data of GSE in this case? I looked for questions that have the same kind of data, but I couldn’t find any, what do you recommend @ATpoint ?

You cannot use GSE data here. It is for microarray as I said above. Download raw data (fastq, this is NGS not array), and then align them. If you are new to this then simply use the FPKM that are provided.

I want to do both from raw data and FPKM @ATpoint, how and where to get the raw fastq data?

On the bottom of the GEO page, you see the "Relations" category, which points to NCBI BioProject PRJNA551397. There, under "project data", you see "SRA experiments". Click on the number 28, which takes you to the sequence read archive page for the samples. Now click on Send to run selector

The recommended way to download these runs from the SRA is to use the SRA toolkit, specifically the program fasterq-dump from it. See here for the How-to. Basically, just run fasterq-dump <RUN ACCESSION>. E.g.: fasterq-dump SRR9606476 for the first run of you experiment. Download the SRA toolkit here

@user432797 You typically download the fastq files either from NCBI or ENA. Check sra-explorer.info which can search and provide download links for datasets. I also wrote a (today a bit deprecated given sra-explorer exists) tutorial that might be interesting, see https://www.biostars.org/p/325010/ It is not recommended imho to directly download files via fasterq-dump, that is often unstable in my experience. The tutorial covers recommended approaches, or simply use the download links from-sra.explorer.info to get fastq right away (my preferred approach).


text:
Discard sequences contain Ns
I have Fasta file contain multiple sequences, some of those sequences contain N's or n's, I'm trying to discard all the sequences and related header that contain Ns or ns in the file.
I have tried multiple python code but return an error.
import sys
from Bio import SeqIO
handle = open(sys.argv[1], "rU")
filtered = [record for record in SeqIO.parse(handle, "fasta") if record.seq.count('N') == 0]
output_handle = open("N_removed.fasta", "w")
SeqIO.write(filtered, output_handle, "fasta")
output_handle.close()
handle.close()

Error:
NsScript.py:3: DeprecationWarning: 'U' mode is deprecated
  handle = open(sys.argv[1], "rU")

and
import sys
from Bio import SeqIO
for record in SeqIO.parse(sys.argv[1], "fasta"):
    if record.seq.count('N') == 0:
        print(record.format("fasta")

Error:
/RemoveNs.py", line 6

                                        ^
SyntaxError: unexpected EOF while parsing

Need to see what the errors are really 'sys.argv' means the file is from the command line

@Michael I have added the error kindly have a look.

The 1st code snippet is not an error, it's a deprecation warning. You can remove that warning by dropping the "U" and just specifying "r" as the mode. (The "U" is for universal newlines, which is now the default behaviour since python 3, and should be specified by a separate newline arg anyway).
The 2nd code snippet is a trivial error because you're missing a closing parantheses ) on the final line.

Second code ... of course and ah yes!

If you're not married to Python, here's a one liner in bioawk:
bioawk -c fastx '{ if($seq !~ /[Nn]/) {print ">"$name; print $seq}}' old_file.fa > new_file.fa

With Visual Studio Code you can delete them using a regular expression. Replace ^>.*\n[^>]*N[^>]* with an empty string
You can also do that in a similar way using the sed Linux command. You can use regular expressions in Python too.

so using command line by applying "" Sed ^>.*\n[^>]N[^>] MyfastaFile.fa  "" will remove them?

I don't know the exact syntax for the sed command... I think it's something like  sed 's/^>.*\n[^>]N[^>]//g' yourfile.fa > yourfilewithoutn.fa

thank you, but this command produces a similar header number of the original file nothing change I assume

That will work in VS Code... I did not try with sed but if you write it with the right syntax it will work too.

There was a  typo in the regex expression due to not marking it as code. It's ^>.*\n[^>]*N[^>]*. The sed command also is wrong for that reason.


text:
Add tags to the read pairs in a bam file
I have a bam file from paired-end sequencing, in which each R1 has the RX tag assigned (via UMItools group) denoting a Unique Molecular Identifier (UMI).
I want to call molecular consensus reads using fgbio but fgbio GroupReadsByUmi gives me an error, because it expects both R1 and R2 to have the RX tag set.
I would like a tool that will go through my bam file and assign each R2 the same RX tag as the R1.

For something that specific you'll probably need to write something.

I would also ask this of the umitools developers: https://github.com/CGATOxford/UMI-tools/issues

The best thing would be to figure out how to do this in UMItools. I don't know that tool at all, but it seems that it is supposed to be able to deal with / tag both reads in paired data. e.g. did you use --paired flag, etc.
I would suggest investigating that route first before you look at the next part of my answer.

Per Devon Ryan's comment suggestion I would suggest a workflow similar to this:

sort bam file by query name (samtools sort -n) to get the R1s immediately following the R2s
Write a short script/command that

reads through the bam file and for each line
looks up the tag in the R1, prints R1 line
reads the R2, substitutes the RX tag in, prints R2 line

pipe that through samtools view and possibly samtools sort back to a bam file of the same type that you started with.

My pysam is a little rusty but something like this might work
import pysam

bam = pysam.AlignmentFile("bam", "rb")
# print header
print(bam.header)

tag = None
for read in bam:
    r2 = False
    if tag:
        read.assign_umi_tag() # assign to read- not sure where/how
        r2 = True
    else:
        tag = read.lookup_umi_tag()  # i don't know where this tag is in bam format
    print(read)
    if r2:
        tag = None

I haven't tested that in any way, caveat emptor.
I am not personally familiar with UMItools and how it assigns tags so can't suggest anything more specific, sorry.

Thanks. My pysam is non-existent, but I can probably do it in Rsamtools.


text:
Fastq dataset after clustering
I'm doing a project for DNA archival storage systems
Is there an open-source dataset used as ground truth for clustering?
E.g
A file containing the ground truth strand sequences, which have been synthesized:
Dataset 1
st1: ATGCGTGC
st2: GCTGATTC
st3: AATGCGTT

And another file containing the sequenced strands, where for each original strand we read multiple 'erroneous' strands:
Dataset 2
st1 cluster:
ATTCTTGC
TTGCGTGC
ATGCGTGT

st2 cluster:
GCTGATTC
GATGAATC
GCAGATTA

st3 cluster:
AATGCGTT
ATTGCTTT
AATATGTT

Where the clusters can then be used for a consensus alignment algorithm.
What I'm aiming to do is to use an algorithm to learn a mapping from dataset1 to dataset2, so that for an arbitrary strand (stx) I can generate a corresponding 'erroneous' cluster (stx cluster)
The two datasets should come from a real synthesis and sequencing pipeline.

I could be misunderstanding the question, but shouldn't just about any barcoded multiplexed sequencing run have this property (in the index read)? How many sequences do you want in your dataset 1? how long should they be? what are the sequencing technologies you are targeting? error models vary dramatically from technology to technology.


text:
UTR features in gbk files
I have a gbk files for a viral genome that I have used biopython to add features of the 3' and 5' UTRs. However when looking in a genome viewer (IGV) I don't see these features. I am not unsure if I encoded the UTR features into the gbk file correctly.
They look like this for example:
5'UTR           39..104
3'UTR           961..2956
5'UTR           1098..1126
3'UTR           2781..2956
5'UTR           3138..3193
3'UTR           6594..7013
5'UTR           15611..15770
3'UTR           16996..17054
5'UTR           complement(17856..17960)
3'UTR           complement(17164..17243)
5'UTR           complement(19557..21099)

Whole gbk file here
I encoded them in the following way:
   seqID = some-gene_name
   Start = SeqFeature.ExactPosition(some_number)
   End = SeqFeature.ExactPosition(some_number)
   seqStrand = 1 #or -1
   FeatureLoc = SeqFeature.FeatureLocation(Start, End)
   Feature = SeqFeature.SeqFeature(FeatureLoc, type = "5'UTR", strand = seqStrand, id =seqID)

I have also tried converting the gbk into a gff file and the line concerning the UTR's look like:
GQ994935.1  feature 5'UTR   39  104 .   +   .
GQ994935.1  feature 3'UTR   961 2956    .   +   .
GQ994935.1  feature 5'UTR   1098    1126    .   +   .
GQ994935.1  feature 3'UTR   2781    2956    .   +   .
GQ994935.1  feature 5'UTR   3138    3193    .   +   .
GQ994935.1  feature 3'UTR   6594    7013    .   +   .
GQ994935.1  feature 5'UTR   15611   15770   .   +   .
GQ994935.1  feature 3'UTR   16996   17054   .   +   .
GQ994935.1  feature 5'UTR   17856   17960   .   -   .
GQ994935.1  feature 3'UTR   17164   17243   .   -   .
GQ994935.1  feature 5'UTR   19557   21099   .   -   .
GQ994935.1  feature 3'UTR   18576   18590   .   -   .
GQ994935.1  feature 5'UTR   21051   21099   .   -   .
GQ994935.1  feature 3'UTR   18576   20039   .   -   .

This is what I see in IGV:

THis is just the regular CDS features with no UTRs
Am I doing this incorrectly?

Summary from comments
The gff3 file has been inspected (presented above), however it remains uncertain how the UTR regions should be represented. The intention is for the UTR regions represented as a features the same as a CDS.
The information copied and pasted is an excerpt from the gbk file, with supporting information from GFF files as requested (below). The formatted described is the Biopython output.
Further information is here: http://pastebin.com/8jm2LEk2

What file are you loading into IGV? how are you expecting to see the UTRs? What do you actually see instead? Please [edit] your question and clarify.

@terdon I have looked at the gff3 file, I am unsure what to expect as I don't know how utr regions are represented but my wish is to see the urt regions represented as a features the same as a CDS. Instead I just see nothing but the normal CDS features

Please [edit] your question and i) explain that you are loading a gff3 file into IGV; ii) show us a few lines of this file so we know what you are doing and iii) include a screenshot of what you see in IGV; iv) clarify if you are working with genbank files or gff files. What are the first lines in your question? They look like they come from a genbank file, but those lines alone are no known format. Is that what you're using?

@terdon. I have tried with both gbk and gff3. What is pasted is just an excerpt from the gbk file (I will add gff also). The format is just that produced by biopython so I am unsure about it being "no known format" as I don't know what the format should look like. Yes that is what i am using :)

I still don't understand what file you are loading in IGV. Is that gff file what you load? The one that only contains UTRs? Or do you load the first file you show, the incomplete genbank file? Can you please give us an example of the actual file you load into UGVIGV?

@terdon I am have tried viewing the gbk file in IGV which is what I have shown above. the gbk and gff contain all the original CDS features but I have just added in the UTRs. here is the complete gbk file I am using: https://pastebin.com/8jm2LEk2

The UTRs display properly in SnapGene or in https://cuba.genomefoundry.org/plot-sequence-features so the problem must be an IGV setting

@Peter Youa re correct, seems to be a IGV error


text:
Can I eject the flow cell MinION from USB port when doing basecalling?
Can I eject the flow cell MinION from USB port when doing basecalling? Because its taking long and its not possible to wait for its finished, will it affects basecalling?

You can, you will lose the transfer from the MinION. It would probably be better to just stop the run. The computer will continue to basecall untill it is done. All the raw reads are stored locally on the machine so it should be fine.

Yes, it should be fine, especially if the run is finished. MinKNOW should still continue basecalling even after the MinION is disconnected. Unfortunately I can't find good posts in the ONT community about that because search is currently broken. This post might be useful (found by searching my email records), but it's two years old:

Hi, thank you for your post. You can certainly remove the flow cell and wash it without affecting the catch-up mode. Please see the downloads page here for the most up-to-date installer version.
- Monolina Binny, ONT

If you're likely to be doing more long runs in the future (>6h), I'd recommend getting a cheap video card with CUDA capability (at least v10) to save time in the future.


text:
Download complete Phylogenetic Tree as JSON
I would like to integrate the phylogenetic tree into a website. Does anyone know where I can download the complete tree? Best would be a JSON, which looks something like this:
[{
  "level": "domain",
  "name": "eukaryota",
  "commonNames": {
    "en": ["eukaryotes"],
    "de": ["Eukaryoten", "Eukaryonten"]
  },
  "children": [
    {
      "level": "kingdom",
      "name": "plantae",
      "commonNames": {
        "en": ["plants"],
        "de": ["Pflanzen"]
      },
      "children": [...]
    },
    {
      "level": "kingdom",
      "name": "animalia",
      "commonNames": {
        "en": ["animals"],
        "de": ["Tiere"]
      },
      "children": [...]
    },
    ...
  ]
},
{
  "level": "domain",
  "name": "bacteria",
  "commonNames": {
    "en": ["bacteria"],
    "de": ["Bakterien"]
  },
  "children": [...]
}]

Other formats and structures are ok, too. The important thing is that the level of the nodes can be retrieved. In the end I only want to have the main levels Domain, Kingdom, Phylum, Class, Order, Family, Genus and Species in the tree. (https://en.wikipedia.org/wiki/Taxonomy_(biology)#/media/File:Biological_classification_L_Pengo_vflip.svg) So if there are more nodes from other levels like Clade, Tribe or Subfamily, it must be possible to identify and remove them.
Thanks a lot!

NCBI maintains an SQL dump of their taxonomic database
You can see what that contains here
For smaller requests (depending on what you need), you might also consider the ENSEMBL REST API - although I don't know of a way to retrieve the complete tree (shoot them an email): https://rest.ensembl.org/taxonomy/classification/53399?content-type=application/json

It would be worth looking at the Tree of Life web project, it isn't downloadable but its a reference source that can be improved on and updated in your website.
ToL contains trees linking everything from bacteria, archaea through to the complete diversity of eukaryotes, but its a bit old with 2008 from my taxonomic knowledge being the last update for some of the entries, e.g. Plasmodium spp..


text:
Are there computational tools to extract features of DNA sequences?
I am looking for tools to extract features from short DNA sequences. For example, entropy, complexity, GC-content, etc.
I have found the generateFeatures.py script from the PyFeat repo, but is there love a more widely used source code or a standard way to extract features from sequences biopython or similar?
Also, I think that since there are $2N$ encoded bits in a sequence of $N$ nucleotides, we have at most $2N$ independent features that could be extracted.

In addition, I am curios if there are any transformers models for DNA sequences.

Edit:
In addition, you can use DeepHF's feature utils, which can be accessed and consumed as seen here.

Great question! +1

If you want to do this with biopython, the SeqUtils package could be a solution.


text:
Comparing AUCs: Discrimination of same Control from different Test Group - paired or unpaired? miRNA
I want to compare two AUCs using a bootstraping method from pROC package (roc.test).
I have done a classification using Logistic Regression and Serum miRNAs as prediction factors. The Controls are in both cases the same, the test group is different.
1. Cancer (high risk) vs. Control = AUC 1
2. Cancer (low risk) vs. Control = AUC 2

If I want to compare this two AUCs - should I use a paired or unpaired test?
I would tend to using a paired test, because of exact the same Control group (Control groups are correlated).

Not sure this can be answered. What are ROCAUCs? Bootstrapping is usually non-parametric, please confirm. ROC curves are a special case of data integrity used in machine learning, I'm not sure thats what you mean. Bootstrapping might be used for data augmentation, we'd need to know if that was the reason. As a simple answer paired tests are more powerful - so yes. The only time this is questioned is if it doesn't give the answer you are looking for.

@Michael ROCAUC is a AUCs of a ROC Curve - nothing special. You are right with the Bootstrapping method - it is non-parametric.
So actually I want to compare two AUCs. Because in both cases I am using the same control group I am not sure if I need a paired or an unpaired test. In both cases the results are not significant and this was what I was looking for. So I do not really care about the concrete p-value - I just want to use the right test for the described case.

Okay I see. It has to be a paired test without bootstrapping.

@Michael Why without bootstrapping? My sample size is quite low with approx. 35/ group. Hence I do not really can assume nomality. This why I wanted to use Bootstrapping. DeLong would be another possibile method.

Generally its paired. The bootstrapping personally I find concerning. I presume you are using it to augment your data set.
Bootstrapping at best tends to be a controversial statistical operation and could result in the non-significance you are observing. Anyway if its not significant between high and low - thats your answer.
You could of course perform a single ML classifying low, and high cancer risk against control. That would provide a method to understand the associations of miRNA.

Bootstrapping issue, I had a think about this. It depends on how the bootstrapping is done. In a paired sample this (I assume) samples with replacement at each paired position for both ROCs. If this is not the case - the test is junk. However, if this is the case, its actually a good test, because the paired structure of the data is always preserved no matter how many bootstrap samples are performed. So like if a 100 were done and the pairwise statistic (which we don't know what it is - but something like KS test) and its >0.05 for all 100 samples - I don't think there's any doubt thats a robust result.
Again if the two ROCs are bootstrapped independently, in my opinion the test is junk.
Reason The concern of bootstrapping is it causes weird distributions which were not present in the original data. This is only okay if there's lots of experience about interpreting the result. For example, in phylogenetics there's decades of experience interpreting bootstrapping and a consensus about what it means. This often isn't the case for most other areas of statistics.
I note this particular issue is really ramping up the rep score ;0)

I have already done a discrimination between Control vs Tumor (low + high). The AUC is quite good with 0.87. Because my Tumor Group as a whole is strongly shifted towards high risk Cancer, I wanted to also show that there is no difference between the AUC of Cancer (high risk) vs Control and the AUC of Cancer (low risk) vs. Control. This way I can strengthen my results.

I see, I think that bootstrappng will bias towards a 'no-significant' difference result. See what the referees say I guess.

Can you please explain briefly, why you think that a paired test should be used? I just want to understand it. Thank you.

@Mischa answered above. Bootstrapping is an opinionated area of statistics outside phylogenetics.


text:
How to visualize called narrowPeak files in UCSC Genome browser?
I have this file:

and I get this error
I googled and found this solution but it still doesn't work and I get this error

Error File 'GSM2797523_FOXA2_IDR0.02_narrowPeak_try.bed' - line 1 of custom input: variable needs to start with letter '"Peak name=""Somite narrowPeak""" '
Any suggestions?
Many thanks
Found a way!!!
opened it with text editor and not exel as before and added the line:
track type=narrowPeak name="Somite narrowPeak"

Could you please forward the first error and  provide some context to your study for example explaining the data behind each column. Sometimes minimal context questions are answered, but usuaully they are ignored.. The title doesn't make grammatic sense  'visualize called narrowPeak files in UCSC Genome browser?'

This is the error
Error File 'GSM2797523_FOXA2_IDR0.02_narrowPeak.bed' - Error line 1 of custom track: thickStart out of range (chromStart to chromEnd, or 0 if no CDS)

Headers should not contain tabs:
track type=narrowPeak name="Somite narrowPeak" description="Somite narrowPeak"

Ensure that you have NO tabs on that line.

I made sure and now I get this error:
Error File 'GSM2797523_FOXA2_IDR0.02_narrowPeak_try.bed' - Unrecognized format type=narrow line 2 of custom track

The type is narrowPeak, not narrow

track type=narrowPeak Peak name="Somite narrowPeak" description="Somite narrowPeak"
Changed it and still get this error:
Error File 'GSM2797523_FOXA2_IDR0.02_narrowPeak.bed' - Unrecognized format line 1 of file: "track type=narrowPeak Peak name=""Somite narrowPeak"" description=""Somite narrowPeak""" (note: chrom names are case sensitive, e.g.: correct: 'chr1', incorrect: 'Chr1', incorrect: '1')

Where is this extra "Peak" coming from?!? Just use EXACTLY what I wrote with ABSOLUTELY NO CHANGES. If you get an error message then it is because you did NOT follow these instructions.


text:
Trying to create a .bam file without the need for a .sam file
I'm trying to use the code specified in this link to create a .bam file without the need for a .sam file.
Here is the code I'm using:
bwa mem -R "@RG\tID:$SAMPLE\tPL:illumina\tSM:$SAMPLE" \
    $REF_PATH/hg38.fasta \
    mapping/R1.fastq.gz \
    mapping/R2.fastq.gz | \
    samtools sort - | \
    tee mapping/$SAMPLE\.bam | \
    samtools index - mapping/$SAMPLE\.bam.bai \
    2>> log

and I get the following error:
[M::bwa_idx_load_from_disk] read 0 ALT contigs
`[E::main_mem] fail to open file ` '.

The ref path and file locations work when creating and converting a .sam to a .bam file.
Any ideas?

The problem is not in skipping the intermediate .sam format, but already the mapping step.

That link does not contain that code — is there a different link it does come from? (Whatever tutorial it's from would benefit from an update. Writing the BAM file via tee so that you can pipe to samtools index is such a spectacular anti-pattern. Though it's unfortunate that sort is one of the subcommands that does not yet support samtools's new --write-index option.)

Ha, I was wrong, the index is indeed loaded. The reads are not...  The expected log should look like this:
mapping [M::bwa_idx_load_from_disk] read 0 ALT contigs [M::process] read 1006360 sequences (150000291 bp)... ...
Then I would follow @Bastian Schiffthaler's answer to find what's the problem with your reads.

It turns out samtools sort does support --write-index since 1.10 after all. So this can be written more simply as bwa … | samtools sort --write-index -o mapping/$SAMPLE.bam - (plus a little more if you want to use a BAI index rather than ‑‑write-index's default CSI).

This is an error from bwa saying that it can't find your input files:
Here's a checklist:

Check that you have no whitespace after your \ newline breaks. I.e. \ is the very last character
Check that you have no whitespace in any path names
Quote your variables: $REF_PATH -> "$REF_PATH"
Use set -x at the top of your bash script to debug.
In general, read this

Why make this so complicated with the tee command?
I don't think you even can index an unsorted bam, so you have to wait until sorting is done before you index.  Indexing is really short compared to mapping and sorting, you aren't saving much by trying to do both at once.


text:
What is the typical host-to-bug DNA ratio found in nanopore sequencing without amplification?
I'm interested in sequencing a human sputum sample using an ONT MinION without performing any type of whole genome DNA amplification or targeted PCR.
Has anyone found a good reference (or anecdotal evidence) for a range of human to pathogen DNA ratios typically found in sputum samples?
Does this depend on whether the patient is infected by a virus or bacteria, and the species that is infecting the patient?
I am seeking a ballpark estimate. Can I expect to find 10x as much human DNA as pathogen DNA, or 100x-1000x?

My recollection of a conversation I had with Justin O'Grady is that it's about 99% host DNA for unfiltered samples. There are a few ways to reduce that portion prior to sequencing. He has written a paper that mentions a saponin-based host depletion method that is quite effective:
https://doi.org/10.1038/s41587-019-0156-5

Anecdotally, I have been told that it varies from person to person, also might be affected by disease or hygiene in a particular environment. Here is an article dealing with that problem:
https://journals.plos.org/plosone/article/comments?id=10.1371/journal.pone.0076096
Not sure if they have references useful for you in their paper, but they describe a way of enriching microbial or human DNA. It might be helpful.


text:
>My counter is counting genotypic combination occurences more than once, how do I ensure it counts one combination and doesnt count it again?
mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

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

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

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

combination_labels = []
combination_counts = []

for lineage in lineages:

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

    for combination in combinations:

            combination_df = df[list(df)]

            for mutation in combination:
                combination_df = combination_df[combination_df[mutation] == 1]

        #print(combination_df.shape[0])

            combination_labels.append('_'.join(combination))
            combination_counts.append(combination_df.shape[0])

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

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

The input CSV
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0

Output CSV
combination,count,percentage
D614G,87355,90.7084929856806

My above code is used to count all occurences of combinations of spike protein mutations

to determine the genotypes.
However, looking at the values processed into the .csv file, i can see ' D614G' is being counted even when
it is in combination with other combinations.

My question is how do i ensure once something has been counted (a row in the imported .csv file)
it is not counted for a second time?
OR possibly how can I edit this code to prevent counting of singular mutations even when presented within a combination?
Thanks in advance

Nice to see a traditional python question based on pandas. The top of your code is chopped off. Could you inset it please?

Could you also add what a row in your CSV looks like?

Input / output csv?

The input CSV in your Result folder is only two lines per file?

The input file has about 97,000 lines of data - each for a strain of SARSCOV2. The output file has around 54 lines, as genotypes but this should be quite a few less

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

for lineage in lineages:

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

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

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

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

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

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

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

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

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

for lineage in lineages:

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

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

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

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

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

Absolutely amazingly detailed reply, thank you so much! It's worked exactly how I needed it! Not sure how to mark this as answered, as I'm very new to this forum - perhaps a mod could help here?
Thank you again, sincerely!

This might help you with that


text:
Python/Biopython - Replace amino acid residue on MSA with "z" from a list of unaligned positions
I'm trying to programmatically replace a set of amino acid residues on an MSA with a "Z" from a list of unaligned positions. Any ideas on how I could do this?
Input: a list of unaligned positions for each sequence + an MSA
Output: an MSA with those unaligned positions on the aligned MSA replaced with a Z for each sequence.
I have a multiple sequence alignment i'm working with and a list of sites on that alignment that are glycosylated. I was trying to run a script that would take the input of the positions as a list of characters along with the MSA to replace the Nth character with a "Z" to easily visualize the sites. So essentially if I feed it the following sequence
M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K

and positions
2, 5, 10, 11, 14

it should spit back
M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K

Hi @DannyMorales and welcome. This is likely fairly easy to do, but the question is difficult to understand. Any clarity would be helpful.

Hey Michael! Thanks for the reply. I have a multiple sequence alignment i'm working with and a list of sites on that alignment that are glycosylated. I was trying to run a script that would take the input of the positions as a list of characters along with the MSA to replace the Nth character with a "Z" to easily visualize the sites.
So essentially if I feed it the following sequence "M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K" and positions "2, 5, 10, 11, 14" it should spit back "M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K"

Understand, what you need to do is update your question with the information above and then provide the code you have so far. The obvious choice here is biopython AlignIO

Could you please show your code, so that is clear what you tried so far?

The code is
import re

pos = [2, 5, 10, 11, 14]
newseq = []
prot = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        newseq.append('-')
    else:
        if (n + 1 - count) in pos:
            s = re.sub(r'[A-Za-z]', 'Z', s)
        newseq.append(s)
print(''.join(newseq))

output

M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K

Precisely as requested, apologies for the delay.
Happy to explain why it works - its just a couple of Python tricks and a method. It must be run under Python 3 (unlikely to work in Python 2) and probably >3.6.
BTW I normally avoid re for complicated reasons but its easy, hence its here.

As @gaspanic has pointed out re is not only unnecessary it will slow the code (re is very inefficient). It is in fact one of the 'legendary inefficiencies' of Python. Thus removing it ...
# variable declarations omitted
newseq = list()
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        newseq.append('-')
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
        newseq.append(s)
print(''.join(newseq))

Same output.

Another alternative, without the need for re:
seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
pos = [2, 5, 10, 11, 14]

new_seq = []
count = 0
for s in seq:
    if s == "-":
        new_seq.append(s)
    else:
        count += 1
        if count in pos:
            new_seq.append('Z')
        else:
            new_seq.append(s)

print(''.join(new_seq))

Result:
M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K

Thanks @gaspanic. You are right re is very slow and removing is of course the better solution altogether. I came to coding via Perl so I tend to over regex and re is one the downsides of Python.

@M__, to be honest I didn't know re was particularly slow in Python. In any case, I think that if efficiency is an issue, it would also be important to define the final size of the new_seq list beforehand (e.g. [None]*len(seq)), and then replace the Nones as one loops over seq.

Hi @gaspanic its very famous for it. It should be compiled via re.compile but its still a problem.
The inbuilt regex - split replace that sort of thing - is usually used instead not for functionality (it ain't as good at all) but performance.
Predefining the list size - I will look into that in detail, thanks and do some tests and timings. I'm well aware for C and Cython it is important.
If I'm doing re I parallelise it.

... okay I tried predefining the list but didn't noticeably improve in my code. deque is another way in. Would make a good question.

Yes, looks like you're right. See this thread: https://stackoverflow.com/a/63600055/8861903 . I guess I must have inherited that idea from other languages (or I'm just plain wrong).

Thanks @gaspanic those timings were better than mine. I'll take another look later. It could be the data set size, this 'data set' is very small and there could be an overhead.

@M__ I actually tried extending the sequence 100000x, but the difference was still minimal. I even removed the if statements to see if that was the bottleneck, but still not much difference. It appears that the explanation lies in how Python allocates memory to lists :).

Hi @gaspanic as promised I've placed a question here: https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example

Great, thx for taking the time! :)

Hi @gaspanic I've placed a 'final answer' in the thread https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example . This combines an iterator with a list comprehension ... Many thanks for your input.


text:
How to see if a gene is the product of copy number variation?
Copy number variation (CNV) is a phenomenon in which sections of the genome are repeated and the number of repeats in the genome varies between individuals.Copy number variation is a type of structural variation: specifically, it is a type of duplication or deletion event that affects a considerable number of base pairs. I have a set of protein coding genes and I would like to verify them to be the product of CNV. My first idea was to overlap their domains, that are well separated by homopeptide regions but before to do this I would like to know whether more suitable indexes exist to infer this characteristic. Can anyone help?

This is going to be rather difficult to measure, I think, with this kind of data. Where did the data come from? Would it be ok to simply find pairs of genes that are homologous? What do you mean by "product of CNV"? How would a gene be the result of a deletion? Can you give us a specific example of the kind of even that you would like to detect?

The data come from NCBI GenBank. Check for homologous genes and to compare their length could be a good strategy

Do you mean you want to find out if a gene is a result of a gene duplication event in a species' lineage? If yes, then something like OrthoFinder could help. Usually I read CNV in reference to variation within a population genetics study. Please expand your description a bit and clarify what it is you are trying to do with more detail.

I have a set of genes taken from cryptosporidium parvum that show off Low Complexity regions. I would like to have a comparison with orthologs proteins and to see if they are longer than, regardless of the LCRs,  orthologs. I would like to establish a connection between CNV and replication slippage if there is one.


text:
Calculating average mito.percentage for each cluster (seurat)
I have a tricky data set with cells that will have a higher percent of mitochondria genes than "typical" data sets. I would like to look at the mito percentage in each cluster without any filtering (e.g., subset(seurat_object, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 & percent.mito < 10))
FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)

Done!
I would also like to calculate the average percent.mito for each cluster and place the value on the FeaturePlot/TSNE plots. I have no idea how to manipulate a Seurat object to calculate these values but I can probably figure out how to place those values on the graph -- probably with geom_text.
My Hail Mary approach is below:
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = p3, FUN = function(x) mean(p3))
CombinePlots(plots = p3)

My attempts yield the error below; I don't know if my approach is close to correct and therefore, I have spent little time working through the error message.
In mean.default(p3) : argument is not numeric or logical: returning NA
Execution halted

and
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = p3, FUN = function(x) mean(seurat_object$meta.data$seurat_clusters))
CombinePlots(plots = p3)

Error in `[[.Seurat`(x, i, drop = TRUE) :
  Cannot find 'meta.data' in this Seurat object
Calls: lapply ... lapply -> FUN -> mean -> $ -> $.Seurat -> [[ -> [[.Seurat

and
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = seurat_object, FUN = function(x) mean(seurat_object$meta.data$seurat_clusters))
CombinePlots(plots = p3)

Error in as.list.default(X) :
  no method for coercing this S4 class to a vector
Calls: lapply -> lapply -> as.list -> as.list.default```

If anyone can provide some guidance -- you will be inducted into my personal hall of fame. :)
Thank you for your time.
P.S. I am using Seurat 3.1
Update:
md = surate_object@meta.data
md_mean <- ddply(md, .(seurat_clusters), summarize, mean_mito_per=mean(percent.mito))

This was not the approach I was hoping for but I have the values -- which when expanded to the length of the original df, I can insert it into a free slot in the seurat_object.
If anyone has a better approach I would love to hear it.
UPDATE:
Using @PPK's approach
# First the cluster annotation and the tsne embeddings are merged
label.df <- cbind(as.data.frame(seurat_object@meta.data), as.data.frame(seurat_object@reductions$tsne@cell.embeddings))

# using dplyr across to calculate the mean mitochondrial percentage and
# the median tsne values per cluster
label.df <- label.df %>% dplyr::group_by(seurat_clusters) %>% dplyr::summarise(dplyr::across(percent.mito, ~ mean(.), .names = "mean_{.col}"), dplyr::across(contains("tSNE"), ~ median(.)))

# running FeaturePlot and adding the label layer
FeaturePlot(seurat_object, reduction = "tsne", features = "percent.mito") +
  geom_text(data = label.df, aes(x = tSNE_1, y = tSNE_2, label = mean_percent.mito))

Error:
* object '.col' is a function.
ℹ Input `..1` is `dplyr::across(percent.mito, ~mean(.), .names = "mean_{.col}")`.
ℹ The error occured in group 1: seurat_clusters = "0".
Backtrace:
     █
  1. └─global::initial_clustering(BM0H, "BM0H")
  2.   └─`%>%`(...)
  3.     ├─base::withVisible(eval(quote(`_fseq`(`_lhs`)), env, env))
  4.     └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
  5.       └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
  6.         └─`_fseq`(`_lhs`)
  7.           └─magrittr::freduce(value, `_function_list`)
  8.             ├─base::withVisible(function_list[[k]](value))
  9.             └─function_list[[k]](value)
 10.               ├─dplyr::summarise(...)
 11.               └

You can extract the necessary values and add them directly the plot as a second layer using plot + geom_text(). This is very similar to the inner workings of the DimPlot function with label = TRUE but allows you to use anything as label.
# I use dplyr v1.0.2 for all data frame manipulations
# First the cluster annotation and the tsne embeddings are merged
label.df <- cbind(as.data.frame(seurat_object@meta.data), as.data.frame(seurat_object@reductions$tsne@cell.embeddings))

# using dplyr across to calculate the mean mitochondrial percentage and
# the median tsne values per cluster
label.df <- label.df %>% dplyr::group_by(seurat_clusters) %>% dplyr::summarise(dplyr::across(percent.mito, ~ mean(.), .names = "mean_{.col}"), dplyr::across(contains("tSNE"), ~ median(.)))

# running FeaturePlot and adding the label layer
FeaturePlot(seurat_object, reduction = "tsne", features = "percent.mito") +
  geom_text(data = label.df, aes(x = tSNE_1, y = tSNE_2, label = mean_percent.mito))

The labels are positioned at the median position of the cells of a cluster. Typing LabelClusters (without brackets) in the R console will show you the labeling function used by Seurat that contains some additional ways of positioning the labels.
If you want to avoid overlapping labels you can use geom_text_repel from the ggrepel package.

I think this is exactly what I need @PPK. As the code is written, I get the following error: Error: unexpected ',' in: label.df <- label.df %>% group_by(seurat_clusters) %>% summarise(across(percent.mito, ~ mean(.), .names = "mean_{.col}"))," Execution halted. I believe summarise should hold both of the across items??? I have tinkered with it at bit, but have yet to resolve the error.

label.df <- label.df %>% group_by(seurat_clusters) %>% summarise(across(percent.mito, ~ mean(.), .names = "mean_{.col}"), across(contains("    tSNE"), ~ median(.)))  Yields this error:
Error: across() must only be used inside dplyr verbs.  ...                     └─dplyr::across(percent.mito, ~mean(.), .names = "mean_{.col}")` -- any tips?

Sorry for the typo in the answer. As you correctly determined the summarise function was closed too early. I modified the answer accordingly.

The second problem seems to be with some other package you are using overwriting some dplyr function. This happens pretty frequently because the functions have such generic names. I updated the answer with explicit calls to dplyr::group:by, dplyr::summarise and dplyr::across to circumvent this kind of problem. Hope this helps

There is still some issue. I will post the error at the end of my post to preserve the format.

Is your code expecting seurat_clusters to be integers rather than characters?

The code should work for both integers and characters.

I looked for solutions and it may be that the error is caused by the expressions ~ mean(.) and ~ median(.). You could try to replace this with simply mean and median.  For me both work equally well but after reading more documentation the second version may be more correct. If this solves your problem I will update the code in the answer again

Let us continue this discussion in chat.


text:
Get list of urls of GSM data set of a GSE set
I have this GSE dataset ( GSE104279 ) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104279).
I want to make a table with set IDs and ftp urls to use it as a table in galaxy.org
I know that we can use ENA to get specific arrangement:
https://www.ebi.ac.uk/ena/browser/
I tried to get:
SampleID    Group   URL

so I used :
https://www.ebi.ac.uk/ena/browser/view/PRJNA412223
But nothing is showing.
Is there away to get these urls in the arrangement above?

You can use Entrez Direct for this as follows:
esearch -db gds -query 'GSE104279' \
  | esummary \
  | xtract -pattern DocumentSummary \
    -if 'entryType' -equals 'GSM' \
    -def 'NA' -element Accession title summary FTPLink

This will return a table with data similar to this:
GSM2580330      18_Z13_2_d5_Zika        cortical organoids_Zika_5d      ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2580nnn/GSM2580330/
GSM2580329      17_Z13_2_d5_Control     cortical organoids_mock_5d      ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2580nnn/GSM2580329/

To download sequence reads, you should follow links to SRA. Using Entrez Direct you can do this as follows:
esearch -db gds -query 'GSE104279' \
  | elink -target sra \
  | efetch -format runinfo

This will return a comma-delimited table containing SRA identifiers and an FTP path to the SRA data. These FTP paths won't lead you to FASTQ files though. You can pass the SRA run identifiers of the format SRR### (or ERR### or DRR###) to fastq-dump or fasterq-dump tools from the SRAToolkit to download data in FASTQ format.

Thank you @vkkodali, this is so helpful, but I want files with sequence, not just summary, by the way, it was good timing, because I just installed Entrez Direct tool but I did not know how to use it, do you know what sequence should I use for cutadapt for these file, they are human genome –