text:
Setting scale of profile plot from zero for chip seq profile plot
Figure from this paper which shows a profile plot in figure4 first figure.
.This i guess is made using deeptools perhaps. Meanwhile
When i tried to do similar plot  the scale doesn't start from 0. So my question how did they set the scale from 0 or the data is as such in other words their data point starts from 0. In my case it starts above 0. So is there a way to scale or normalize to make it start from 0 as i see in the embo paper or its just the characteristics of the data itself generated.

lol.. your sequencing coverage is super high i guess, such that the average coverage across the genome, even for your control is 10 reads per million? i cannot see a y-axis label

and did you flip the samples? I think you should not bring it to zero.. it is what it is. If your control is a chip input, this is what I would expect. Note the example you showed above, it is a KO.

"and did you flip the samples?" that mystery haunts me since I'm the end guy in the chain it our collaborator who are might have jumbled up the samples i have checked  quite  number of times regarding the sample labeling im sure im running the samples as provided.

"f your control is a chip input, this is what I would expect." Its the  WT ,h2bub chip of control sample

yes for my plot the color are fine ..my doubt was regarding if the scale can be set to zero which i understand is a bad idea

yeah you should not.. maybe use a zscale for the axis. Like convert your coverage to a zscale for each gene, take average z scores across all genes and plot that

you should see a flat line around 0 for your control and a bump for your chip

"maybe use a zscale for the axis. Like convert your coverage to a zscale for each gene, take average z scores across all genes and plot that " it can be done with deeptools or something else required?

i code it from scratch.. you can maybe use https://deeptools.readthedocs.io/en/3.0.0/content/tools/computeMatrix.html and just scale the rows


text:
Picking number of samples for clustering
I have expression data for 100 samples and clinical data for 50 among them. I want to carryout k-means clustering using expression data and then do a survival analysis (Kaplan-Meier analysis) using the labels obtained from k-means.
Can I used the expression data of all 100 samples to carryout clustering or should I use only 50 samples for clustering?

Hi @user98059 I understand your question perfectly at a analytical level, but I think it would be preferable to edit it and provide a biological explanation regarding survival analysis.I'm not sure it will attract many responses. ... anyway. In summary, I would go for 100. It also depends on whether K-means is inside or outside ML methods.

In case there is no systematic difference between the samples with or without the clinical data (i.e. one set is microarray and the other is RNA-seq), I would use all of the samples. And if the cluster labels obtained with "clustering 50 samples" and those obtained with "clustering 100 samples" does not differ too much, then doing either of two would not have much effect anyway. If you would like to compare cluster labels in that way and assuming you would want to test a number of ks, it would be better to compare each combination of ks and number of samples.

Hi @Haci The thing that bothers me is it isn't clear with K-means is within supervised or unsupervised learning. In my experience it can be used for either. These are quite different calculations. At a guess it is unsupervised because in supervised learning all the samples would be used without any question.

source: Stack Exchange

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

sort: votes


text:
Is it possible to filter contaminated reads for raw PacBio sequences (not HiFi reads) before assembly?
De novo genome assembly for non-model organisms face the issue of bacterial contamination. For assembled contigs with mostly bacterial-like sequences (based on BLAST search), the entire contig can be dropped off. But for a long contig with small fragments of bacterial-like sequences, one way to solve it is to break the contigs into smaller pieces and drop the bacterial-like pieces. It would work better if the raw reads could be removed BEFORE assembly. This was commonly done for Illumina reads, but for PacBio long reads (non-HiFi reads) with a higher error rate, is it possible to filter contaminated bacterial sequences based on sequence similarity, like BLAST, Kraken? Of course, the potentially contaminated region should be long enough to drop the entire sequence, for PacBio long reads. Maybe 30%? Any suggestion of the processing details?

I don't think that in principle the issue is any different for PacBio and Illumina; it's just that the process might be noisier with higher error rates. Note that PacBio errors aren't randomly distributed substitutions and probably many (most?) short-ish k-mers (<30nt) will likely be intact. I would suggest doing k-mer filtering of your reads as suggested here or here or here.
Alternately, use your favorite alignment-based technique, if you have a good idea of what your contaminant is. In principle PB reads will be better filtered by this method than Illumina as local alignment will be able to figure out the errors quite well with longer queries.
If you go the route of trying to do alignment to some set of contaminant references, the obvious answer would be to use minimap2 as your aligner as this is exactly the task that it is designed for. You can then use samtools to filter out anything that maps well.
Here is a related question/answer that also addresses the contaminant issue and suggests using Blobology, which I'm not familiar with but may be helpful (granted it also seems to be Illumina-centric). Genome-size estimation using k-mers on PB reads appears to be pretty accurate too, so I think that it's likely ok to work with raw PB k-mers.

Right, the PacBio CLR reads errors are randomly distributed. Using smaller kmer is a great idea. I find the BBsketch is the fastest comparison using kmer analysis, although the documentation is embedded deep in a readme file. BBSketch could build a local NCBI nr sketch and compare reads against it in less than an hour. I am comparing the accuracy of BBSketch and Kraken (micro-organisms) now.  It might be hard to track down the sources for de novo assembly, within unknown bacteria, parasites, and even plants.  Blobology is now Blobtools2, a very nice tool I have been exploring! Thanks!

You can use the kraken database. It will classify your reads based on what organism it matches with. Then can write a custom script to extract the reads you are interested in for downstream analysis.

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


text:
how to train a gene dataset with a nearest shrunken centroid classifier?
I have a data file named "geneexp.csv".
the data contains information about gene expression of three different cell types (CD4 and CD8, CD19) I want to classify cells by performing the nearest shrunken centroid classification of training data in which the threshold is chosen by cross-validation. I split the data (70% train and 30% test).
data = read.csv("geneexp.csv")

splitData <- function(data, trainRate) {
  n <- dim(data)[1]
  idxs <- sample(1:n, floor(trainRate*n))
  train <- data[idxs,]
  test <- data[-idxs,]
  return (list(train = train, test = test))
}

split <- splitData(data, .7)
train <- split$train
test  <- split$test

then with the use of pamr package I tried to buid the following model and plot :
y <- train[[ncol(train)]]
x <- t(train[,-ncol(train)])
mydata <- list(
  x = x,
  y = as.factor(as.factor(y)),
  geneid = as.character(1:nrow(x)),
  genenames = rownames(x)
)

# Training and cross-validating threshold
model <- pamr.train(mydata)
cvmodel <- pamr.cv(model, mydata)
pamr.plotcv(cvmodel)

but I can't make it work. I get the following error:

Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels

I have already transfer y to the factors. Can you help me? How can I fix it?

as.numeric(as.character(y)) this looks incredibly dangerous. What is in y ? can you do table(y)

in y there is 3 classes CD4,CD8,CD19. I got this error: NAs introduced by coercion

table(y)
CD19  CD4  CD8
68   74   68

It's a bit different from usual R, but if you check the help:

data: The input data. A list with components: x- an expression genes
in the rows, samples in the columns)

So in your case, you need to transpose the matrix:
mydata <- list(
  x = t(x),
  y = as.factor(y),
  geneid = as.character(1:nrow(x)),
  genenames = rownames(x)
)

As I don't have your amazing data, I can only use iris below:
pamr.train(list(x=t(iris[,1:4]),y=iris[,5]))

pamr.train(data = list(x = t(iris[, 1:4]), y = iris[, 5]))
   threshold nonzero errors
1   0.000    4       6
2   0.841    4       7
3   1.682    4       10
4   2.523    4       11
5   3.364    4       13
6   4.205    4       18
7   5.046    3       22
8   5.887    3       23


text:
Pairwise alignment using CIGAR String
I want to perform pairwise alignment and get the resulting CIGAR string (Concise Idiosyncratic Gapped Alignment Report). A fuller description being here
The channel can emit insertions, deletions and substitution errors.
For example, if the input was AAGCT
Then the output could be AAAAGGC
Where there were the following errors in the following positions.
1, Insertion of A at first position
2, Insertion of A at second position
3, Insertion of G at the third position
4, Deletion of T at the last position
I want to know how I could obtain the CIGAR string given two sequences. In addition, could I use viterbi algorithm for this?

What's your question?

How to obtain CIGAR string given two sequences

Hi @DukeGlacia could you check the answer? It looks pretty good, and thereafter consider accepting it?

Following the basic rules and preferring the AAG contiguous
we can write:
  01234
  AAGCT     (ref)
AAAAGGC     (query)

So the CIGAR string is 2I3M1I1D that is:

2I insert 2As  at the start:
3M match the 3 bases AAG
1I insert the G
1D delete the final T


text:
chipseeker annotation issue
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
library(clusterProfiler)
library(org.Hs.eg.db)

peaks <- readPeakFile("MACS_folder/pol2_chip/Pol_2pausing/pol2_diff_peak.bed")

print(peaks)

#peakAnno <- annotatePeak("merge/blast.bed", tssRegion=c(-3000, 3000),
#                        TxDb=txdb, annoDb="org.Hs.eg.db")

peakAnno <- annotatePeak(peaks, tssRegion=c(-1000, 1000),
                         TxDb=txdb, annoDb="org.Hs.eg.db")

So this is basically what im doing to annotate the peaks.
Now the issue comes with annotation for example I see TLR4 annotated to this
118082896   118083202   307 *   shp1_R1_peak_679    Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9   117708197   117708979   783 1   7099    ENST00000490685.1   374699  ENSG00000136869 TLR4

Now The peak is 679 as you can see from the IGV plot as well as annotated file. But when i see in IGV browser TLR4 is far away from the peak where i get the annotation.
So now is it an issue with the code im using or due to the inherent nature of the data.
Any help or suggestion would be really appreciated.

I had a similar case just recently.
The tssRegion argument ChIPseeker is used to define what is considered a promoter, which is frequently the question at hand. However, ChIPseeker will not perform any kind of filtering on the annotated peak file.
In the results there is a column annotation that contains information what kind of genomic region the peak is in. If it is within your tssRegion window it will be labeled as a promoter.
In your case the annotation is 1st Exon = Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9. This tells you that the peak is not part of the promoter region you defined and located within the first exon of the indicated transcript. The peak is really far away from TLR4 (374699 bp) but I think all other sequences nearby are pseudogenes and this is why TLR4 is selected as the closest gene.
You can get an overview over the annotations using plotAnnoPie(peakAnno)
tl;dr: you need to filter the results of the ChIPseeker annotation to get the genomic regions you want.

Glad this is a reproducible problem and now its more confusing for some peak the genes assigned are perfectly fine matches  to the coordinate given but for others its would be in the intro or some exon but would show as annotation being promoter. " ChIPseeker annotation to get the genomic regions you want."  how shall i do that ? simple distanetoTSS filter?

you can either directly use the distance to TSS to or define the promoter region via tssRegion and then filter all peaks where the annotation contains "Promoter" i.e. grepl("Promoter", annotated.df$annotation)

i did annotate again with homer it seems pretty accurate as compared to chipseeker ran homer with default parameters it works fine ..need to find why the same can;t be seen with chipseeker


text:
What is the best way to map/align my reads on a given genome?
I am frequently using a ballgown package for my rnaseq analysis, but recently I have had a new task to have my reads mapped on two different genomes to understand the level of alignment between the two, therefore I wonder what is the best tool to do this? I have single-end reads.
What I need is to have a percentage of uniquely mapped fragments, mapped to multiple loci and, say, unmapped. Can bowtie do the job? Is there a standard explanation procedure explained elsewhere? Thanks.

This is a pretty straightforward analysis, assuming that you want to measure something like "how many reads align to genome 1, and how many reads align to genome 2?".

The bowtie2 manual has a section on this.
There are other tutorials for bowtie2.
bwa is also well-suited to this task.

You should be able to then use samtools to compute statistics on your alignments.


text:
MUMmer plot error. Line 884. What to do?
My third year comp sci friend got me this far, I hope you guys can help me with the last inch! For the record, I have practically no coding experience. I am using the 3.9.4 alpha version of MUMmer.
Here is the manual I have been following: http://mummer.sourceforge.net/examples/
I am at part 2.2.2 Running the Mummer plot, and I type in the final line to generate the plot.
I type:
mummerplot -x "[0,8666124]" -y "[0,9025688]" -postscript -p mummer mummer.mums
And the error is:
Can't use 'defined(%hash)' (Maybe you should omit the defined ()?) at /usr/bin/mummerplot line 884

Line at 884:
if ( defined (%$rref) ) {                                                                                #-- skip reference sequence or set atts from hash
        if ( !exists ($rref->{$idR}) ) { next; }
        else { ($refoff, $reflen, $refdir) = @{$rref->{$idR}}; }

Its perl

Hi could you please present your code inline, e.g. copy and paste? It is hard to read as an image. This looks like Perl, could you confirm this?

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

I had the same error and I solved by removing defined and writing only if (%$rref) in the mummerplot code. You have to do this for all the lines that you have it. There should be six lines to correct in this way .
You can find an explanation of why here, that tells that the function defined:

Returns a Boolean value telling whether EXPR has a value other than the undefined value undef. If EXPR is not present, $_ is checked.

And then:

Use of defined on aggregates (hashes and arrays) is no longer supported. It used to report whether memory for that aggregate had ever been allocated. You should instead use a simple test for size:
if (%a_hash)   { print "has hash members\n"   }

@KamilSJaron thank you but to be honest the code is written in C and I am not really expert on this. Anyway, I tried to add more details that is what has helped to me .

This is totally enough for people to get the idea that it's just an outdated part of the code, good job in figuring it out!


text:
Snakemake wildcard constraints reject too much
I have the following two rules:
rule extract:
    input: "data/{dir}/{file}.gz"
    output: "data/{dir}/{file,^.*(?<!\.gz)$}"
    conda: "config/conda-extract-env.yml"
    shell: "cd 'data/{wildcards.dir}'; gunzip -k {wildcards.file}.gz"

rule extract_dot:
    input: "data/{dir}/wtdbg2.3.dot.gz"
    output: "data/{dir}/wtdbg2.3.dot"
    conda: "config/conda-extract-env.yml"
    shell: "cd 'data/{wildcards.dir}'; gunzip -k wtdbg2.3.dot.gz"

The second one is actually a hack because the first does not work on files named wtdbg2.3.dot.
I have also tried giving the constraint using the wildcard_constraints member of the rule, and prefixed the string with an r or not.
Whatever I do, as soon as I uncomment the extract_dot rule, snakemake complains about the missing input file data/C.elegans/wtdbg2.3.dot in another rule. Note that there is a rule with output data/C.elegans/wtdbg2.3.dot.gz, and when changing this rule to output data/C.elegans/wtdbg2.3.dot instead, everything works as well.
So I conclude that the problem is in the regex. I have tested the regex on regex101.com, and it accepts wtdbg2.3.dot. But snakemake seems not to. Does anyone know if there is anything special about the snakemake regexes? Should they not be normal python regexes?

Does anyone know if there is anything special about the snakemake regexes? Should they not be normal python regexes?

No, nothing special.  The Snakemake documentation just links to the Python re module documentation.
I think this is something strange in how the anchors ^ and $ behave when matching wildcards, and not anything else in your particular regular expression.  Check out the pattern-matching with this dummy Snakefile:
rule all:
    input: ["data/subdir/example", "data/example", "example"]

rule anchored_example_of_subdir:
    """Can't add either ^ or $ or this will fail"""
    output: "data/{subdir,subdir}/example"

rule anchored_example_in_subdir:
    """Can't add ^ or this wil fail"""
    output: "data/{file,example$}"

rule anchored_example:
    """This works with both ^ and $"""
output: "{file,^example$}"

(There are no inputs or actions so it's just a matter of what filenames match what rules, and the regular expressions do nothing but match exact strings.)  This "works" as is (snakemake -nrp shows what you expect) but if you add in anchors as described in the docstrings above the upper two rules no longer match the expected filenames.  It's as if (and maybe is?) under the hood Snakemake is trying to apply the expression to the entire path.
Can you just drop your ^ and $?  The snakemake documentation implies that it's using the equivalent of re.match rather than re.search already, with wording like "Here we could restrict the wildcard dataset to consist of digits only using \d+ as the corresponding regular expression" so I don't think it's necessary anyway.

Edit: I think I found the issue in the snakemake source.  It munges your constraints together into one larger regular expression for the entire file path, so ^ and $ get embedded inside the resulting expression.  This looks like something that should at least be clarified in the documentation.

Thank you for your answer! I will report this as an issue then.


text:
shiny code with seurat object
I am trying to create a shiny board for my scRNA seq data. I have saved my seurat objects adn I use following code as below:
library(shiny)
library(shinythemes)
library(Seurat)
library(uwot)

set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")

ui <- fluidPage(

    # Application title
    titlePanel("scRNAseq data browser"),

    fluidRow(
        column(4,

               sidebarPanel(
                   radioButtons("dataset", label = h3("Dataset"),
                                choices = list("sample1" = "set_1", "sample2" = "set_2"),
                                selected = "set_1"),
                   helpText("Enter gene names in CAPS.  Gene names must be exact."),
                   textInput("gene2", label = "Gene Name", value = "GAPDH")

               )
        ),
        column(8,
               column(6,
                      mainPanel(
                          plotOutput("dimPlot1", width = "150%"),
                          plotOutput("genePlot2", width = "150%")
                      )
               )
        )
    )
)

server <- function(input, output, session) {

    datasetInput <- reactive({
        df <- (input$dataset)
            })

    output$dimPlot1 <- renderPlot({
        DimPlot(datasetInput(), reduction = "umap")
    })

    output$genePlot2 <- renderPlot({
        FeaturePlot(datasetInput(), features = (input$gene2), reduction = "umap")
    })

}

shinyApp(ui, server)

That error that I get on the section where image needs to be displayed is "Subscript out of bound".
Following dataset can be downloaded and used as an example curl -Lo pbmc_10k_v3.rds https://www.dropbox.com/s/3f3p5nxrn5b3y4y/pbmc_10k_v3.rds?dl=1
Can anyone help in fixing this shiny script so that it works.
Many thanks.

Not easy to reproduce your error but here is my guess. This is a string, you can check it:
input$dataset
class(input$dataset)

This is a simple solution. If you want to call out the dataset depending on input, use a list:
set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")
data_list = list(set_1=set_1,set_2=set_2)

Then in the server :
datasetInput <- reactive({
        df <- data_list[[input$dataset]]
            })

It worked, thank you:)

@StupidWolf if you can give a look at this https://bioinformatics.stackexchange.com/questions/21043/shiny-app-for-single-cell-cluster-export-issue


text:
Troubles implementing LIMMA for paired samples (before/after treatment) comparaisons
i try to implement LIMMA for paired samples in order to compare gene expression before/after treatment. But...i'm not confident in my results since every single gene expression seems to be significantly associeted with the treatment!
Here is my methodology:
mat<-t(df)
where df is a dataframe with:
-1000 columns for different gene expression
-+ 1 column 'id' (10 different id, repeated 2 times)
-+ 1 column 'ttt' (0=before, 1=after treatment)
-20 rows of observations (10 patients observed twice)
So once I get my matrix:
subject_id<-df$id A vector with 10 different id repeated 2 times
ttt<-as.factor(df$ttt) A vector with the treatment status (before/after)
Here is my design matrix:
design<-model.matrix(~ 0+ttt)
colnames(design) <- levels(ttt)
As asked in limma() description i specify the correlations:
corfit <- duplicateCorrelation(mat,design,block=subject_id)
Fiting of the model:
fit<-lmFit(mat,design,block =subject_id,correlation = corfit$consensus)
fit2<-eBayes(fit)
Getting adjusted pvalues for every gene comparison before/after treatment:
topTable(fit2, adjust.method='BH',number=1000)
But every single gene expression seems to be relevant, so i'm very doubtful about my methodology
Any idea?
Many thanks!

what do you actually have, RNAseq or microarray data? I suspect you have RNAseq data and you need to normalize it

Hello, thank you for your response. I have transcriptomic responses, with logarithmic transformation (one value per gene expression).

What u have above should be ok. Just note that for the last step u are testing all coefficients to be zero. Not sure if that is what u want

But if you specify treatment and groups correctly. It should be ok. Without more info, hard to tell


text:
Finding a gene of similar length
I would like to find a few genes that are similar in length to a gene of interest, and by that I mean the total length of a canonical transcript including introns. This is proving to be surprisingly hard using web based lookups.
The gene of interest is 14.5Kb long and I wanted a table of genes and their canonical lengths to compare and contrast to in order to pick a few to have a look at in some vcfs I have.
I have a low bioinformatics ability but happy to give things a go!
All the best
I have tried:

first download a file of all gene files from ensembl which looks like this:

Then in R:
x <- read.delim("all_genes_bed.txt")
x$length <- (x$Gene.end..bp.-x$Gene.start..bp.)
x <- subset(x, x$length >14000 & x$length <15000)

Then I just had a look in the RStudio browser.
That seemed to work!

I would simply download a GTF annotation files, filter for type=gene, subtract end from start coordinate to get the length, and from this make a simple lookup table, listing gene name and length. This you can sort by length and then check genes that intersect the length of your gene. Do you need code suggestions?

@ATpoint many thanks for taking the time to answer. Would I get this from ensembl? I think once I have the file I could do that in R. Not sure how to get the whole gene table though

@tactolimus Yes, you could get it from Ensembl. How about trying something, see how far you get and if you hit the wall then just edit it to the toplevel post and we can have a look to debug it.

@ATpoint I think I solved it, let me know if this seems reasonble

Yes, I think this is reasonable.


text:
How do I change the identity of a sample from spatial 'Visium' data preprocessing in Seurat v3?
I'm wanting to create a merged object in Seurat using 2 10x Visium 'slices'. However, when I create an object, Seurat assigns an identity "SeuratProject" to the objects (by default I'm guessing; see attached plot). This will be problematic when I attempt to merge the objects as there will be unintended duplicates of the objects.
Here's what I'm attempting to use to create the 10x Spatial Object:
brain1 <- Load10X_Spatial(data.dir = "~/xxx", filename = "filtered_feature_bc_matrix.h5", assay = "Spatial", slice = "slice1", filter.matrix = TRUE, to.upper = FALSE)

I'd need the ability to assign unique ident to each of these (and any subsequent) slices/data.
Thank you!!

You can re-assign orig.ident after creating the object. For example:
brain1 <-Load10X_Spatial(data.dir = "~/xxx", filename = "filtered_feature_bc_matrix.h5", assay = "Spatial", slice = "slice1", filter.matrix = TRUE, to.upper = FALSE)
brain1$orig.ident <- "Visium"

Didn't work for me. I don't know what I did wrong.

The active.ident is still "SeuratProject", so also add Idents(brain1) <- "Visium".  There should be a "project" option in Load_10X_Spatial to change the "SeuratProject" default


text:
Mapping protein refseq to Gene ID
I have a protein refseq (eg, NP_000029). How can I get the corresponding gene ID (ag, APC) from NCBI using an R package?

Hi @ElizaGanguly Is  efetch not appropriate here?

You can do this with the reutils package, which provides an API to NCBI's E-utilities. Here's an example for your specific question:
install.packages("reutils")
library(reutils)

# Get universal identifier
uid <- esearch("NP_000029", db = "gene")
# Fetch summary
sm <- esummary(uid, db = "gene")
# Extract specific tag
sm$xmlValue("//Name")


text:
Analyzing and visualizing metagenomic data
I am starting out in metagenomics,
Aim I wish to understand the whole analysis from the raw reads to interpretation.
Background My current pipeline is used for:

analyzing human gut microbiota cleans reads,
filters out human genome, and
maps the reads to known microbial species using IGC and mOTUs.

Question I am looking for basic tools to visualize the output, in order to get some biological and statistical insights (Linux and python are preferable environments).
For example, an obvious step would be to display the mOTUs analysis results on the phylogenetic tree, using the full phylogeny provided with this package. Whilst I could implement it myself, I was wondering if anyone else had previous experience?
Update
One package that I found useful is MOMR (through in R).


text:
GC normalization ATAC seq data
Is it really required to do GC normalization ATAC seq data ?
One of the paper where they have ATAC seq data did the normaisation of GC bias after peak calling. When i looked for the library and the reason why GC bias i found this
Could it be true that genes with higher GC content are higher expressed? It has been suggested that genes that are either extremely high or extremely low expressed are under some form of selection leading to “extreme” GC content. What CQN does, is making the effect of GC content comparable across samples, and we show in [1] that this leads to improved inference. It also flattens the effect of GC content on gene expression, but we believe this is better than having the effect of GC content depend on the sample.
https://www.bioconductor.org/packages/release/bioc/vignettes/cqn/inst/doc/cqn.pdf

Now that scenario is for RNA seq would it be logical to do the same for ATAC seq?

It is not required if you ask me. GC content is the same across samples per genomic region so it would perhaps make sense if you want to compare regions within the same sample, but then other biases such as mappability would need correction as well, not sure if the bother is necessary. I personally do not do it. It is not established from what I know and I do not see the advantage. You can simply plot GC content vs peak counts to explore yourself if you like.

" I personally do not do it. It is not established from what I know and I do not see the advantage. " i have the doubt but not the knowledge  if i should do or not ! ..i came across paper where they did that but not sure why ..but now im clear how to proceed thank you and @Devon Ryan for the answer

You will always find papers doing this or that. It comes down to whether the paper shows conclusively whether what they did is beneficial or necessary, and to my knowledge this has not been done for ATAC-seq towards GC normalization, whereas for RNA-seq it has been shown (see Alpine package at Bioconductor, Mike Love the developer/first author).

The only reason to normalize for GC content in RNA-seq is if it differs notably between samples/groups. If that's not the case and you aren't trying to compare genes withing samples then you have no reason to try to account for GC content. The same goes for ATAC-seq, though there the danger with trying to normalize for GC content is that you then mask changes in global accessibility.
In short, normalize for things that should cause problems with your analysis and nothing else.

Now my confusion is cleared


text:
Summing up one column entries for each sample
I have a dataframe with 27 samples and for each of them I have different number of values
I want to sum up the values for each sample, leaving each sample with a single value
> dput(d[1:10,c(1,8)])
structure(list(SampleID = c("LP6008460-DNA_A04", "LP6008460-DNA_A04",
"s30", "s30", "s67", "s15", "s67", "s67", "s67", "s30"), X.1 = c(8.17407e-05,
1.40067e-06, 3.86017e-05, 9.50247e-05, 7.08333e-07, 3.19133e-05,
3.1602e-05, 5.00513e-05, 3.29613e-05, 1.003e-06)), row.names = c(NA,
10L), class = "data.frame")
>

I started to do that manually but that was a pain and inaccurate

I have answered this question, but it seems to be more suitable for stack overflow since it's just a programming question without any bioinformatics in. And it's also best to at least pretend you've made an effort to answer the question before, since this is a very common question and already has many answers.

It's worth learning how to use dplyr for this sort of thing:
dat %>%
    dplyr::group_by(SampleID) %>%
    dplyr::summarise(sample_mean = mean(X.1))

# A tibble: 4 x 2
  SampleID          sample_mean
  <chr>                   <dbl>
1 LP6008460-DNA_A04   0.0000416
2 s15                 0.0000319
3 s30                 0.0000449
4 s67                 0.0000288

Sorry for each sample should be sum not mean

I think you can make that tiny change yourself, @Exhausted

Thank you @RamRS I replaced mean with sum myself but did not make sense

What do you mean by "did not make sense"?


text:
How to deal with the mismatches between gene names obtained from different sources?
For most of the time, I rely on gene ids to combine different datasets. However, in some instances, I have to combine datasets based on gene names. Then, if I don't know the source of gene names in the dataset, I get to this issue of choosing a source of gene names, be it Ensembl, HGNC etc for human genes. I wonder if this is a common issue and if there is an reliable method out there to deal with this issue.
To demonstrate the mismatch between different sources, I compared gene names for all human genes. I obtained them from 4 different sources as listed below, using BioMart (pybiomart) :
+-----------------+-----------------------+--------------------------------------------+
|     source      |     attribute_name    |                display_name                |
+-----------------+-----------------------+--------------------------------------------+
| HGNC            |  hgnc_symbol          |  HGNC symbol                               |
| NCBI            |  entrezgene_accession |  NCBI gene (formerly Entrezgene) accession |
| Uniprot         |  uniprot_gn_symbol    |  UniProtKB Gene Name symbol                |
| Ensembl (maybe) |  external_gene_name   |  Gene name                                 |
+-----------------+-----------------------+--------------------------------------------+

Upon this comparison, I found several things that are clearly apparent.
1. Genes names of protein coding genes are best matched across different sources.
I saw that protein coding genes have the best matching (left, measured in terms of Jaccard index) across different sources, with majority of genes having a single unique names (shown on right).

However, there isn't a good enough matching in the case of not protein coding genes. Here, HGNC and Ensembl have the best match. (I don't expect Uniprot gene names to match because they are of course only for protein coding genes.) Remarkably most of the genes have 2 unique ids (shown on right).

2. Gene names from some databases match with each other.
Comparison of all genes shows that some pairs of the sources do not have a good match e.g. Ensembl and Uniprot, with many genes having 2 unique gene names(!).

I saw similar pattern for genes on chromosomes (autosomes,X,Y) and on the scaffolds.

3. Mitochondrial gene names do not match at all(!).
Mitochondrial genes clearly have different names in different databases. None of the genes have a single unique gene names (!).

How to deal with such a mismatch between different sources?
Should I prefer one particular source or is there a way to make use of the synonymous gene names from different sources?

This is one hell of a challenge and unfortunately, there's no "right" way. All of these databases are evolving at different rates and the 80+% overlap you see is the best one can hope for. Most LINC/LOC/MT/-DT/-AS1 gene symbols will differ between sources and even if you were to use the most reliable identifier (which IMO is the ENSG identifier from EnsEMBL), you'll still need to manually sort through stuff and use coordinates and google searches to "unify" these gene names/symbols/identifiers.

It seems that Consensus CDS (CCDS) project (https://www.ncbi.nlm.nih.gov/projects/CCDS/CcdsBrowse.cgi) aims to unify gene names/ids of protein coding genes between Ensembl, NCBI and HGNC.

Good catch. But then, CCDS focuses on protein coding regions only, which are already stable (maybe CCDS has a role to play in that already). The other genes cannot be easily cross-mapped.

Here is a post that might interest you: https://www.biostars.org/p/478165/


text:
Query on htseq count
This question has also been asked on Biostars
I am trying to run htseq-count for carrying out rna-seq analysis for solanum tuberosum and i used the following command:
htseq-count --format bam --order pos -s no -a 10 -t exon -i gene_name --idattr gene_id SO_8612_L1.bam GCF_000226075.1_SolTub_3.0_genomic.gff > L1_htseq_count.tsv
and im getting an error message saying:
[E::idx_find_and_load] Could not retrieve index file for 'SO_8612_L1.bam' Error processing GFF file (line 12 of file GCF_000226075.1_SolTub_3.0_genomic.gff): Feature exon-XM_015312074.1-1 does not contain a 'gene_id' attribute [Exception type: ValueError, raised in features.py:326]

How to resolve this? and also can gff be used instead of gtf file. what is the standard gtf file used for Solanum tuberosum?


text:
How to extract certain genes including non-coding regions?
I would like to look at the non-coding regions of some genes. I have around 1000 full genome assemblies and I was able to extract the nucleotide sequences for certain proteins with Roary. However, we are also interested in comparing non-coding regions. Is there a tool that can do that or a tool that tells me what region in the full assembly is relevant? The assembly and extracted genes are present as fasta and fa.ali (alignment) files. The genomes are from bacteria.

LoL, sure, but I don't what all genes, only blaZ for example.

I found Piggy. https://academic.oup.com/gigascience/article/7/4/giy015/4919733

Get a GTF/bed file of the gene coordinates and use bedtools perhaps

Can you give more information about what you did with Roary?  I'm not familiar with it, but the docs say it takes GFF3 as input, so could your input files already have what you need?

Roary is a pan-genome tool. I processed full genome assemblies first with Prokka which generates GFF3 files and then with Roary to get a gene presence/absence map. Yes, Roary uses GFF3 files, which probably are indexed. That is a good point. Is there a tool though that can extract certain genes from GFF3 files? Including non-coding part.

I don't know enough about what existing tools are out there that might do it for you without extra coding, but the format is really just tab-separated text, so at the very least you can read them as spreadsheets and see what the feature column looks like.  For example in some GFF3 files I have from NCBI I see things like region, gene, mRNA, exon, and CDS, each with associated start and end positions.  So you could take the region for, say, "gene," to get both coding and non-coding parts.  (Depends on what's in your files, though!  Maybe share a few lines as an example.)


text:
Cromwell: can't find file ID
I am using Cromwell to test a CWL (common workflow language) implementation of the GATK somatic CNV (copy number variation) pipeline. I am currently getting the following error (full log at the end):

Cannot find a tool or workflow with ID file:///var/folders/cm/0j97xqy577z2dlfht90hv_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl in file /var/folders/cm/0j97xqy577z2dlfht90hv_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl

The ID's are the same except for the first starts with "file://". I am not very familiar with Cromwell so I'm not sure how to proceed in resolving this error. Any help would be appreciated.
Thank you
Full log:

[2020-12-03 15:20:24,82] [info] Running with database db.url = jdbc:hsqldb:mem:28b7f7a3-d261-42ee-b941-e329cdcf46ba;shutdown=false;hsqldb.tx=mvcc
[2020-12-03 15:20:29,97] [info] Running migration RenameWorkflowOptionsInMetadata with a read batch size of 100000 and a write batch size of 100000
[2020-12-03 15:20:29,98] [info] [RenameWorkflowOptionsInMetadata] 100%
[2020-12-03 15:20:30,03] [info] Running with database db.url = jdbc:hsqldb:mem:ae8f4f13-d319-4439-8681-b684e13ca6ac;shutdown=false;hsqldb.tx=mvcc
[2020-12-03 15:20:30,23] [info] Slf4jLogger started
[2020-12-03 15:20:30,25] [info] Pre Processing Workflow...
[2020-12-03 15:20:30,27] [info] Pre Processing Inputs...
[2020-12-03 15:20:30,46] [info] Workflow heartbeat configuration:
{
"cromwellId" : "cromid-d8084bc",
"heartbeatInterval" : "2 minutes",
"ttl" : "10 minutes",
"writeBatchSize" : 10000,
"writeThreshold" : 10000
}
[2020-12-03 15:20:30,48] [info] Metadata summary refreshing every 2 seconds.
[2020-12-03 15:20:30,50] [info] KvWriteActor configured to flush with batch size 200 and process rate 5 seconds.
[2020-12-03 15:20:30,50] [info] WriteMetadataActor configured to flush with batch size 200 and process rate 5 seconds.
[2020-12-03 15:20:30,50] [info] CallCacheWriteActor configured to flush with batch size 100 and process rate 3 seconds.
[2020-12-03 15:20:30,78] [info] JobExecutionTokenDispenser - Distribution rate: 50 per 1 seconds.
[2020-12-03 15:20:30,80] [info] SingleWorkflowRunnerActor: Submitting workflow
[2020-12-03 15:20:30,82] [info] cwl (Unspecified version) workflow 837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a submitted
[2020-12-03 15:20:30,82] [info] SingleWorkflowRunnerActor: Workflow submitted [38;5;2m837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a[0m
[2020-12-03 15:20:30,82] [info] 1 new workflows fetched
[2020-12-03 15:20:30,82] [info] WorkflowManagerActor Starting workflow [38;5;2m837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a[0m
[2020-12-03 15:20:30,83] [[38;5;220mwarn[0m] SingleWorkflowRunnerActor: received unexpected message: Done in state RunningSwraData
[2020-12-03 15:20:30,84] [[38;5;220mwarn[0m] Couldn't find a suitable DSN, defaulting to a Noop one.
[2020-12-03 15:20:30,84] [info] Using noop to send events.
[2020-12-03 15:20:30,85] [info] WorkflowManagerActor Successfully started WorkflowActor-837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a
[2020-12-03 15:20:30,85] [info] Retrieved 1 workflows from the WorkflowStoreActor
[2020-12-03 15:20:30,85] [info] WorkflowStoreHeartbeatWriteActor configured to flush with batch size 10000 and process rate 2 minutes.
[2020-12-03 15:20:30,86] [info] MaterializeWorkflowDescriptorActor [[38;5;2m837d7e7d[0m]: Parsing workflow as CWL v1.0
[2020-12-03 15:20:30,91] [info] Pre-Processing /var/folders/cm/0j97xqy577z2dlfht90hv_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl
[2020-12-03 15:20:32,24] [[38;5;1merror[0m] WorkflowManagerActor Workflow 837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a failed (during MaterializingWorkflowDescriptorState): cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor$$anon$1: Workflow input processing failed:
Cannot find a tool or workflow with ID file:///var/folders/cm/0j97xqy577z2dlfht90hv_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl in file /var/folders/cm/0j97xqy577z2dlfht90hv_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor.cromwell$engine$workflow$lifecycle$materialization$MaterializeWorkflowDescriptorActor$$workflowInitializationFailed(MaterializeWorkflowDescriptorActor.scala:200)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor$$anonfun$2.applyOrElse(MaterializeWorkflowDescriptorActor.scala:170)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor$$anonfun$2.applyOrElse(MaterializeWorkflowDescriptorActor.scala:165)
at scala.runtime.AbstractPartialFunction.apply(AbstractPartialFunction.scala:34)
at akka.actor.FSM.processEvent(FSM.scala:665)
at akka.actor.FSM.processEvent$(FSM.scala:662)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor.akka$actor$LoggingFSM$$super$processEvent(MaterializeWorkflowDescriptorActor.scala:123)
at akka.actor.LoggingFSM.processEvent(FSM.scala:801)
at akka.actor.LoggingFSM.processEvent$(FSM.scala:783)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor.processEvent(MaterializeWorkflowDescriptorActor.scala:123)
at akka.actor.FSM.akka$actor$FSM$$processMsg(FSM.scala:659)
at akka.actor.FSM$$anonfun$receive$1.applyOrElse(FSM.scala:653)
at akka.actor.Actor.aroundReceive(Actor.scala:514)
at akka.actor.Actor.aroundReceive$(Actor.scala:512)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor.aroundReceive(MaterializeWorkflowDescriptorActor.scala:123)
at akka.actor.ActorCell.receiveMessage(ActorCell.scala:527)
at akka.actor.ActorCell.invoke(ActorCell.scala:496)
at akka.dispatch.Mailbox.processMailbox(Mailbox.scala:257)
at akka.dispatch.Mailbox.run(Mailbox.scala:224)
at akka.dispatch.Mailbox.exec(Mailbox.scala:234)
at akka.dispatch.forkjoin.ForkJoinTask.doExec(ForkJoinTask.java:260)
at akka.dispatch.forkjoin.ForkJoinPool$WorkQueue.runTask(ForkJoinPool.java:1339)
at akka.dispatch.forkjoin.ForkJoinPool.runWorker(ForkJoinPool.java:1979)
at akka.dispatch.forkjoin.ForkJoinWorkerThread.run(ForkJoinWorkerThread.java:107)

[2020-12-03 15:20:32,24] [info] WorkflowManagerActor WorkflowActor-837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a is in a terminal state: WorkflowFailedState
[2020-12-03 15:20:37,26] [info] SingleWorkflowRunnerActor workflow finished with status 'Failed'.
[2020-12-03 15:20:40,51] [info] Workflow polling stopped
[2020-12-03 15:20:40,52] [info] Shutting down WorkflowStoreActor - Timeout = 5 seconds
[2020-12-03 15:20:40,52] [info] Shutting down WorkflowLogCopyRouter - Timeout = 5 seconds
[2020-12-03 15:20:40,52] [info] Shutting down JobExecutionTokenDispenser - Timeout = 5 seconds
[2020-12-03 15:20:40,52] [info] Aborting all running workflows.
[2020-12-03 15:20:40,52] [info] JobExecutionTokenDispenser stopped
[2020-12-03 15:20:40,52] [info] WorkflowStoreActor stopped
[2020-12-03 15:20:40,53] [info] WorkflowLogCopyRouter stopped
[2020-12-03 15:20:40,53] [info] Shutting down WorkflowManagerActor - Timeout = 3600 seconds
[2020-12-03 15:20:40,53] [info] WorkflowManagerActor All workflows finished
[2020-12-03 15:20:40,53] [info] WorkflowManagerActor stopped
[2020-12-03 15:20:40,53] [info] Connection pools shut down
[2020-12-03 15:20:40,53] [info] Shutting down SubWorkflowStoreActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] Shutting down JobStoreActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] Shutting down CallCacheWriteActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] SubWorkflowStoreActor stopped
[2020-12-03 15:20:40,53] [info] Shutting down ServiceRegistryActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] Shutting down DockerHashActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] JobStoreActor stopped
[2020-12-03 15:20:40,53] [info] Shutting down IoProxy - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] CallCacheWriteActor Shutting down: 0 queued messages to process
[2020-12-03 15:20:40,53] [info] WriteMetadataActor Shutting down: 0 queued messages to process
[2020-12-03 15:20:40,53] [info] KvWriteActor Shutting down: 0 queued messages to process
[2020-12-03 15:20:40,53] [info] CallCacheWriteActor stopped
[2020-12-03 15:20:40,53] [info] DockerHashActor stopped
[2020-12-03 15:20:40,53] [info] IoProxy stopped
[2020-12-03 15:20:40,53] [info] ServiceRegistryActor stopped
[2020-12-03 15:20:40,54] [info] Database closed
[2020-12-03 15:20:40,54] [info] Stream materializer shut down
Workflow 837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a transitioned to state Failed
[2020-12-03 15:20:40,57] [info] Automatic shutdown of the async connection
[2020-12-03 15:20:40,57] [info] Gracefully shutdown sentry threads.
[2020-12-03 15:20:40,57] [info] Shutdown finished.

I updated to the newest version of Cromwell and the error went away!

The question was by the same person, so I'm guessing they're providing an update.

Hi @Stephanie Hoyt, if this answered your question, could you please mark it as the accepted answer?


text:
Metagenome simulation from a concatenated FASTA file
I am trying to simulate metagenome sequencing. So I will start with a file with a lot of concatenated genomes. From there, I would like to randomly extract 10,000 sequences of length 200bps.
I got help subtracting a random 200 sequence from a fasta file like this:
>JQ086372.1 Enterobacteria phage HK446, complete genome
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTCGCAGAAAGAGATTTTTCACCTTAAAGGCTTCGGATTCACCGGGCTTGTAGGCCTGTCAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

using this:
tail -n +2 file.fasta | tr -d '\n' > newfile
n=$(stat -c "%s" newfile)
r=$(shuf -i1-"$((n-200+1))" -n1)
newfile tail -c+"$r" | head -c200

And as an output I got: GTTAAGGGCGCCGTTCTGCTTGTCTTTGCTGACATGTTTGAACACCGGACGGCACAAAGCGAAGTGCAGCTTTATGAGAATGCAGCCGCCGAACGCATGATGTTCATTCATCGCAACTGGCGCGGTAAATCTGAACCTGAGGAGGGCTCCTGATGGAACCTGGACGATTCAGGCACAGGGTAAAAATTCTCACCTTCACG
But now I want to extract instead of 1 fragment of 200 characters; 10,0000 random fragments of 200 characters using the strategy above. For this, I am looking to exclude all those lines that begin with ">" (p.e. >JQ086372.1 Enterobacteria phage HK446, complete genome).The input file looks like this:
>JQ086372.1 Enterobacteria phage HK446, complete genome
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTCGCAGAAAGAGATTTTTCACCTTAAAGGCTTCGGATTCACCGGGCTTGTAGGCCTGTCAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

>NC_010391.1 Salmonella phage Fels-1, complete genome
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGAGCGTTTGAGTGCCCCATTTGTGAAGCAATAAAGTTCGGGTTCGCGCCAGCGGCAAGCGCCCAGCATGCACCGATTTTTTTAGTTTTAAGGTGTTGACATTAGGTATGTCGGGACTGGTATGCTTTCCTGTGTCGCAGCCCGGCGCGTCTCAATGCAGATTCCCATATCCTGTTCAGGGAAGGTGCAGCATAATGTGCTTTGGATTTACCTGAGCGATCGGTAATTTGCGGATTGAAGACAAAGGTAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

>NC_007577.1 Prochlorococcus marinus str. MIT 9312, complete genome
AATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTACAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC

So for my new output, I will end up with a file with 10,000 subsets of 200 characters that would look like this:
>1
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTTCAGTCTGCCAATATGCGTTCCGACAGCAGATAAAC

>2
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGCAGATA

>3
AATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAAT

Makes sense? Any help?
Thanks in advance

This is a great question for discussion. @gringer  is 100% correct, the question is which and why

It might be better to formulate your question, or append the summation of the goals of your study to match them against a simulator

Given that you want to simulate short-read sequencing from reference genomes, there are read simulator programs that already exist. For example, see this answer.

Its a fantastic discussion


text:
Genes associated with digestion in single cell RNASeq
I am reading Saunders et al. (2018).
In the methods' subsection "ICA based analysis and clustering" (independent component analysis), a set of independent components from ICA is labelled as "Artifact" with the following explanation:

Artifact (n = 315), in which either: (a) cell scores showed obvious replicate-specific biases, or (b) the top loading genes largely intersected with those observed to vary strongly with digest time (Figures 3 and S3), or contained many heat-shock proteins;

Aside from a possible typo (Figure S3 instead of S2) I couldn't find a list of genes that vary with digest time. The list of genes they found in few technical components is shown in Fig. S2 (D-H). The legend of the figure indicates:

Example of digest-related ICs from FC co-analysis of under- and over-digested cells. Left, cell loadings color-coded by digest condition. Right, gene-loading plots for the 10 top and bottom loading genes. IC 15 (cluster 2, Syt6+ deep-layer glutamatergic) and IC 68 (cluster 1, Nptxr+ upper-layer glutamatergic) show contributions from similar genes. Genes related to ATP synthesis and Calmodulins (eg., Cox6a1, Ndufa4, Calm2) load on the under-digested cells, while another set of genes including the nuclear-enriched Meg3 and Malat1 load highly onto the over-digested condition.

However, a full list of genes associated with digestion is not provided. I searched for other articles reporting such genes but I was unable to find any. Is there such a resource or well-known characterization of genes associated with technical signal in scRNA-Seq?


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

from PIL import Image

    def segmentImage():
        segment_image = instance_segmentation()

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

Cell example image:

Example segmentation result:

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

Wow, what an honor that you personally respond. Thank you very much. I actually already found this article and I wanted to avoid labelling the cells by hand. But I guess it pays off in the long run hopefully. I see a lot of potential in your work for Bioinformatics in general and especially for the project i'm working on now. If done right I'm sure PixelLib can help save us a lot of work :)

Thank you......


text:
What does the 'suppressed' folder mean on the NCBI ftp genomes website?
I wish to understand the term suppressed in folder nomenclature on the NCBI ftp genomes website?
Example: here ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/Oscillatoriales_cyanobacterium_C42_A2020_001/all_assembly_versions/
Normally, I'd see folders with names like representative or latest assembly, which are pretty straightforward. But what does suppressed refer to in this case?
Thank you for your time!

Here you have the explanation:

Suppressed genomes: Several hundred RefSeq bacterial genomes that did
not pass assembly or annotation quality validation have been
suppressed(RefSeq genome quality criteria). An FTP report file listing
the most recent suppressions will be provided soon in the RefSeq FTP
site.


text:
Parsing PLINK recombination genetic map
I have a Plink genetic map with 4 columns: Chromosome, rs-snp, Map(cM), Position(bp) and approximately 8,000 lines find the recombination rate  [Rate(cM/Mb)] from a HapMapII genetic map of 4 columns: Chromosome, Position(bp), Rate(cM/Mb), Map(cM) ~ 47,000 lines ... some chromosomes will be 6 times larger than this; not all Plink entries will be found in HapMapII, but 99% will.
I need to parse PLINK output and spreadsheet VLookUP is too slow.
Basically I have to find one missing field matching two given fields in a much larger file. I have tried variations of grep on Linux and FindStr in windows but they have not worked.
Thank you!
Here is the Plink map
Chromosome, rs-snp, Map(cM), Position(bp)
19  rs2312724   0.00187 266034
19  rs1020382   0.002681    267039
19  rs12459906  0.010664    276245
19  rs10409452  0.068867    288246
19  rs3813154   0.092406    292082
19  rs895344    0.099705    292901
... 8300 more rows

Here is the HapMapII map
Chromosome  Position(bp)    Rate(cM/Mb) Map(cM)
chr19   253938  0.154179    0.000000
chr19   256859  0.154528    0.000450
chr19   259772  0.154530    0.000900
chr19   260970  0.154697    0.001086
chr19   261033  0.154888    0.001095
chr19   266034  0.347963    0.001870  <-
chr19   266544  1.256622    0.002047
chr19   266888  1.333840    0.002480
chr19   267039  1.542843    0.002681
chr19   267213  2.171166    0.002950
chr19   267614  1.108918    0.003820
chr19   269380  0.711702    0.005779
chr19   276245  1.791471    0.010664
chr19   277717  2.004476    0.013301
chr19   277776  3.133699    0.013420
chr19   279645  4.135431    0.019277
chr19   280299  4.235347    0.021981
chr19   281360  6.164471    0.026475
chr19   287129  6.113884    0.062038
chr19   288246  6.135026    0.068867
chr19   288320  6.136509    0.069321
chr19   292082  8.817796    0.092406
chr19   292583  9.059132    0.096824
chr19   292901  9.062522    0.099705
.... over 46,000 more lines

I need the Rate(cM/Mb) field imported into a version of the Plink map.
Position(bp) 266034 ; and Map(cM) 0.00187 exist in both files; but Rate(cM/Mb) 0.347963 does not. I need to import Rate(cM/Mb) from the larger file into a new version of the smaller file. The columns can be resorted. I don't need the chromosome column. I will add a final column with mutation rate as well.

What you want to do is a join on two columns, these answers show you how to do it. Is that enough for you or I can write a full-length answer later?


text:
Genomic relationship matrix explanation
I am seeking the definition of the genetic distance matrices used in genomic prediction.

What is the difference between a Rodger distance and kinship matrix?
what people mean when they talk about genomic relationship matrix as in this paper:
https://www.g3journal.org/content/9/5/1545#ref-30 , viz.

then we obtained the Cholesky decomposition of the genomic
relationship matrix

You may want to take a look at this question (and the answers) in the biology community: https://biology.stackexchange.com/q/35330/59521


text:
Calling variants with DeepVariant on targeted NGS sequencing (custom library)
I am seeking advice regarding DeepVariant analysis.
To avoid false positives I'm using several variant callers and then the resulting common set will be considered as TP variants.
One of the callers is DeepVariant (I'm using it by the first time). I used DeepVariant-inception_v3-1.0.0+data-wes_standard deep learning model (checkpoint) as a argument in call_variants function. However I don't have typical WES data (whole exome sequencing), but custom/target NGS pointed to plenty potential regulatory sequences i.e. promoters, TFBS and also pseudogenes and other non-coding stuff.
I have doubts whether my approach is correct, hence many variants may be discarded.
Should I use WGS (whole genome sequencing) model instead of WES model regarding training set preparation?
Or maybe this tool isn't a good choice regarding such custom experiment?
Best Regards, Adam


text:
ATAC seq peak annotation : Which peak to choose for downstream analysis
After doing peak call I'm annotating the peaks using chipseeker tool, which I want to take further downstream analysis. Now unlike RNA seq where i have single gene and their respective expression. But here the peaks for a single gene multiple peaks. As an example have filtered promoter peaks.
seqnames  start    end width strand     V4 annotation geneChr geneStart geneEnd geneLength geneStrand geneId
1     chr1 826797 828101  1305      *  Peak1   Promoter       1    826832  852225      25394          1 643837
2     chr1 869647 870324   678      *  Peak4   Promoter       1    868240  870201       1962          2 284593
3     chr1 876462 877140   679      *  Peak5   Promoter       1    868071  876903       8833          2 284593
4     chr1 877219 878105   887      *  Peak6   Promoter       1    874529  877234       2706          2 284593
5     chr1 923459 924318   860      * Peak14   Promoter       1    923928  939291      15364          1 148398
6     chr1 924379 926062  1684      * Peak15   Promoter       1    925150  935793      10644          1 148398
       transcriptId distanceToTSS         ENSEMBL    SYMBOL                                    GENENAME
1 ENST00000623808.3             0 ENSG00000228794 LINC01128 long intergenic non-protein coding RNA 1128
2 ENST00000432963.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
3 ENST00000446136.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
4 ENST00000635557.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
5 ENST00000420190.6             0 ENSG00000187634    SAMD11    sterile alpha motif domain containing 11
6 ENST00000437963.5             0 ENSG00000187634    SAMD11    sterile alpha motif domain containing 11

Now if we see first 6 rows we have two genes with three different peaks. When i load and see in the igv browser I do see the there are subtle differences in the accessibility signals.
How to address this issue ,in other words which peak to consider or take all the peak and consider their accessibility as one.
Any suggestion or help would be really appreciated.

I'll consider all peaks around promoter of particular gene as one peak. Moreover I'll interpret this as an open chromatin region = active gene. Why? Why there's no one wide peak which represents open chromatin in promoter of particular gene? Why is divided into several narrower peaks? There might be some transcription factors or other proteins binded to DNA, hence you see peak (open, free DNA) - break (valley, TF binded) - peak(again free DNA where Tn5 cut DNA and incorporate adapters) in IGV.
In my study I compared ATAC-seq and H3k27Ac. In such case usually you may notice that peaks from chip and atac never literally overlap. Because ATAC represents open chromatin (DNA without any binded proteins) and chip histone modification (DNA is actually wound on histone proteins), hence we have something like this figure C, first we see high ATAC peak - open chromatin, then H3k27ac peak - lysine acetylated = gene activated. In other words ATAC-seq won't show you peaks around histones, binded polymerase, or transcription factors.

excellent answer to clarify my doubt ..this your paper i suppose will go through it https://bioinformatics.stackexchange.com/questions/15041/atac-seq-peak-annotation-which-peak-to-choose-for-downstream-analysis/15042#15042 to find out enhancer " ATAC-seq reads should map to known histone modifications associated with these features." this needs to be done then only one can comment on the enhancer feature if i understand

You can see from your table that the different peaks correspond to different promoters for different transcripts, for example, FAM41C has ENST00000432963.1, ENST00000446136.1, ENST00000635557.1. And the fact that they are all within the TSS (see distanceToTSS) suggest these peaks are actually falling onto actually TSS.
The underlying issue is are these TSS or do will we mis classify the peaks that are inside genes or other non-coding elements as TSS peaks. Nothing to do with TF, histones etc.
I cannot see the number of reads in those peaks, but if one of them is really strong, and the rest almost non-existent it would make sense.
How can we improve this, I think it depends on the aim of your analysis. Since you are looking at TSS, I guess it makes sense to be able to link it to something that is actually expressed. You can do:

To use an updated annotation (i see that some of the transcripts come from old hg37 annotation)

you have prior knowledge which transcripts might be expressed from RNA seq data, I would use that to reduce the number of transcripts

excellent insight

"To use an updated annotation " i used chipseker hg38 annotation

Ok.. yeah, then check whether all those have the same number of reads.. and try to subset on transcripts that are expressed. Then you can sum up all the reads, or use only the top expressed transcript


text:
How to generate a consensus sequence from a multi-reference BAM file?
I am trying to generate a consensus sequence from a BAM file that was generated by mapping reads to a reference FASTA containing multiple sequences.
Usually, I generate consensus sequences from BAM files using samtools and bcftools:
samtools mpileup -vf reference.fasta filename.sorted.bam | bcftools call -m -O z - > filename.vcf.gz

bcftools index filename.vcf.gz

bcftools consensus -f reference.fasta filename.vcf.gz > filename.consensus.fasta

However, this doesn't work when the reference FASTA contains more than 1 sequence.
Is there a way to generate a consensus sequence from a multi-reference BAM file?
It would also be awesome if the method filled in gaps in coverage with Ns or Xs in the generated consensus sequence but this isn't essential.
Thank you for your help!

This post on biostars has several suggestions.
THe one that seems most straightforward, starting where you are, is to use vcftools:
cat ref.fa | vcf-consensus file.vcf.gz > out.fa

Based on that I think that you should also be able to solve this problem using alfred and GATK.
If all of these suffer the same problem of needing a single sequence (which I very much doubt that they do), you could split the bam / ref files into multiple files, run whatever consensus workflow, and then combine them again. But that doesn't seem like it should be necessary.
See also this post and this post.


text:
Adding matrix as a confounder in EnrichmentBrowser
Brain data, particualrly human brain data have very heterogeneus cell types, so it's important to normalise by cell type/add proportoins of different cell types to the formula when performing differential expression. We have some bulk RNA-seq brain data from disease samples and wildtypes. I have used the package BrainInABlender to get a matrix of the proportions of each cell type in each sample

Cell type
Sample_1
Sample 2

Neuron
0.5
0.7

Microglia
0.5
0.3

The rows are the types of cells, the columns are the samples. I wanted to know what's the best way to add this as a confounding variable in Enrichment Browser (or if that's not possible in DESeq2).
Edit: Current way I've tried to correct for this is in DESeq2. I have transposed the matrix of cell type proportoins then set DESeq's formula:
~ Condition + Condition:Neuron:Microglia:etc

This has drastically reduced the number of genes found to be differentially expressed but I don't know if it's the right way to do it.
Edit 2: When I tried removing the interaction terms, replacing them with +
~ Condition + Neuron + Microglia etc

As suggested on the Bioconductor support section the number of genes went up by an order of magnitude from just using condition which is the opposite of what should have happened. This persisted whether condition was at the start or end of my formula.

What have you tried? DESeq2 works with count data not decimal numbers. What is your biological question/goal?

I have no idea what is an enrichment browser. If you are referring to RNAseq analysis,  in DESeq2, when you set up the data with colData = ... just add the cell type neuron / microglia as a numeric column?

@llrs Sorry for not being clear eough. I have a seperate matrix of count data. I also have this matrix which has the proportoins of each cell type in each sample. My goal was to best correct for differences in proportoins of cell types in each sample.

@StupidWolf Enrichment Browser is an R package that allows for running several differential expression methods easily. I have added it to the coldata and seperated with : in my formula (Condition:Neuron:Microglia:etc but have no idea if that's right.

@Sethzard this doesn't solve the last question. Correct for differences in proportion of cell types for what? To see what cell type express?

@llrs Different types of cells express genes at different levels. Brain samples vary significantly in the proportions of cell types from one sample to the next. A sample containing a higher proportion of one cell type may report changes associated with differences in cell types rather than true differential expression unless samples are normalised or changes in those proportoins are corrected for. We have brain samples from a control group and a disease group and want to work out the truly differential expressed genes. I've edited the post to try to clarify too.


text:
10x Single-cell RNA-Seq - Updating gene symbol with latest ensembl release - Ensembl IDs missing
I am trying to update the gene symbols provided by the 10x human precompiled reference (2020-A), using the mapping Ensembl ID to gene symbol provided in release 102.
However, I found out that 56 of the Ensembl IDs (which I presumed to be stable) are not present in the latest release. See the attached image.

See, for example, gene OCLM (ENSG00000262180) which seems to have a retired identifier:

Is there a way to automatically filter or correct these identifiers, or should it be done manually? What is the most correct approach to take in this case (e.g., ignore them)?
Update
All but one are associated with deprecated identifiers, e.g., ENSG00000262180
The only one with 3 currently mapped identifiers is ENSG00000183791.
Is it correct at this point to filter them out because associated with deprecated identifiers? What about the one with currently mapped identifiers?


text:
Why GAG protein HIV-1 (virus) is in Staphylococcus sp. SS21 and in Klebsiella quasipneumoniae (Bacteria)?
I am one ignorant...
Why GAG protein HIV-1 (virus) is in Staphylococcus sp. SS21 and in Klebsiella quasipneumoniae (Bacteria)?
These are the links in GenBank:
Staphylococcus sp. SS21
LOCUS       WP_193573829             503 aa            linear   BCT 29-OCT-2020
DEFINITION  hypothetical protein [Staphylococcus sp. SS21].

HIV-1
LOCUS       KU749412                9188 bp    RNA     linear   VRL 08-DEC-2016
DEFINITION  HIV-1 isolate DEMC14PK009 from Pakistan, partial genome.

protein align:
mgarasllsggqldrwerirlrpggkkkyrlkhivwasrelerfavnpgllessegcrkileqlqpslqtgseelrslyntvavlycvqqridvkdtkealekieeeqnkskkkaqqtaadkgnnspvsqnfpivqnlqgqmvhqpisprtlnawvkvveekafspevipmfsalsegatpqdlntmlntvgghqaamqmlketineeaaewdrlhpvhagpvapgqmreprgsdiagttstlqeqigwmtnnpaipvgeiykrwiilglnkivrmysptsildirqgpkepfrdyvdrfyktlraeqasqdvknwmtetllvqnanpdcktilkalgpaatleemmtacqgvggpshkarilaeamsqvtqptnilmqrgnfrtprktvkcfncgkeghiakncraprkkgcwkcgkeghqmkdcserqanflgkiwpshkgrpgnflqnrpeptappappeesfrfgegtttpaqkqepidkdlyplaslkslfgndpssq
|||||| | |  || |||||||||||| | ||| ||||||| ||| ||||||  |||  |  |||| |||| |||||| |||| ||||   | |||||||| ||||||||   || |  |        | | ||||||||||||| |||||||||||| ||||||||||||| ||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||||||||||||| ||| | | ||  ||||||||||||||||||||||| |||||||||||||||||||| |||||||| |||||||| ||||||||||||||| |||| |||||||||||||||| |||| ||||||| |      ||| ||     | |||||||||| | ||||||||||||||||||||||| ||||||||||||||||||||||| |||||||||   ||||| | || || |||  ||   || ||||||| || ||
MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNKCQQKAKQ--AEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQTT--NAIMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPA---ESFRF-EETT-PAPKQEQKDK--EPLTSLKSLFGSDPLSQ

The links in GenBank of Klebsiella quasipneumoniae :
LOCUS       WP_199660587             458 aa            linear   BCT 06-JAN-2021
DEFINITION  hypothetical protein, partial [Klebsiella quasipneumoniae].
ACCESSION   WP_199660587
VERSION     WP_199660587.1

protein align:
-----------------------grkkyk-lkhlvwasrelerfalnpgllesaegcqqimeqlqsalgtgseelkslfntvatlwcvhqrieikdtkealdkveevkskqkaqqaqrtqqaaaatgsssqnypivqnaqgqmvhqsmsprtlnawvkvveekafspevipmftalsegatpqdlnmmlnivgghqaamqmlkdtineeaaewdrthpvhagpvapgqireprgsdiagttstlqeqigwmtnnppipvgeiykrwivlglnkivrmyspvsildirqgpkepfrdyvdrffktlraeqatqevknwmtesllvqnanpdcktilrglgpgatleemmtacqgvggpghkarvlaeamsqaqqsnimmqrgnfrg-qrtikcfncgkeghiarnckaprkrgcwkcgkeghqmkdcnerqanflgriwpshkgrpgnfpqsrpeptappaesf-geevtpsprqe-----pvd--kel-----ypp
                       | ||   ||||||||||| ||||||||||  ||| ||  ||| || || ||| ||||||||| |||  |  ||||||||  ||   |   | |   |  ||| |  | ||||||| |||||||  ||||||||||| |||||||||||||||||||||||||| ||| |||||||||||||||||||||||| |||||||||||| ||||||||||||||| ||| | | ||| |||||||||| |||||||||||| ||||||||||||||||||||||||||||||||||||||| ||||||||||||||| |||||||||||||||||||||||||||||||||     ||||| || |  ||||||||||||| |||| |||| |||||||||||||||||||||||| |||||||||||| ||||||||||||||  || || | ||     |    | |
MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNK-CQQKAK--QAEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQ-TTNAIMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPAESFRFEETTPAPKQEQKDKEPLTSLKSLFGSDPLSQ

The links in GenBank of hypothetical protein HMV63_24500 [Escherichia coli] and protein align with HIV-1 GAG:
LOCUS       HAJ5901352               246 aa            linear   BCT 11-MAY-2020
DEFINITION  TPA: hypothetical protein HMV63_24500 [Escherichia coli].
ACCESSION   HAJ5901352
VERSION     HAJ5901352.1

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------mtsnppipvgdiykrwiilglnkivrmyspvsildirqgpkepfrdyvdrffktlraeqatqdvknwmtetllvqnanpdcktilralgpgasleemmtacqgvggpshkarvlaeamsqannanimmqksnfkgprrnvkcfncgkeghiarncraprkrgcwkcgkeghqmkdcterqanflgkiwpshkgrpgnflqsrpeptappaesfkfeettptprqeprdrepltslkslfgsdplsq
                                                                                | ||| ||| ||||||||||||||||||| |||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||| |||||||||||||| ||||||||||||  || |||| |||||| |  |||||||||||||||||||| ||||||||||||||| |||||||||||||||||||||||||||||||||||| |||||| | ||| | |||||||||||||||||
MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNKCQQKAKQAEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQTTNA-IMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPAESFRFEETTPAPKQEQKDKEPLTSLKSLFGSDPLSQ

I remember Blast find result in plant bacteria or virus too

for me could be one lab escape...
it seems the page is for me this
  https://science.sciencemag.org/content/266/5193/1981

there is a pay wall so I can not read that...
The title is
Crystal structure of the catalytic domain
of HIV-1 integrase:
similarity to other polynucleotidyl transferases

i know nothing of argument, but it seems someone unknow for unknow reason infect one bacteria of Klebsiella pneumonia with the virus HIV-1. Someone know what happen if one common bacteria infect with the HIV-1 is diffused? what happen if there is something activate virus replication in the bacteria when the bacteria is in some human?
the genbank file would be:
LOCUS WP_142766724 415 aa linear BCT 07-JAN-2020
DEFINITION DDE-type integrase/transposase/recombinase, partial
 [Klebsiella pneumoniae].

The general speech
We, You don't have the right (or the intelligence) to change
one iota of any DNA or RNA. Genetic material, if we still
want to live a few years, must always be read-only,
never re-written from us because consequences are unknows.

Eukaryotic viruses cannot infect bacteria. The cell wall is thick. The membrane has different protein (entrypoints) and then even if you transformed a bacterium with a viral genome on a plasmid it would not express because of promoters and very different codon bias. And... if you made a codon-optimised viral proteome under bacterial promoters and transformed it the assembled viral particles would not be able to get out.

for me not deep enought

The interesting numbers of an alignment are coverage (91%) and identity (82.68%) of API70200.1 (the protein encoded in KU749412) against WP_199660587.
There are two cases:

technical error: samples were contaminated, which happens not infrequently.
biology did something cool: horizontal gene transfer

The former is likely to be present when you have a single sequence where this gene is present and it is not from a complete genome.
The latter, unless it's an ancient transfer, would likely result in flanking regions and genes.
So blasting the gene to the Klebsiella quasipneumoniae proteome, shows there are three hits (two poor). There is no strain name in list, whereas there are three complete genomes of three strains. So the sequence WP_199660587 is a one off and not from a complete genome. Therefore, the clues strongly point towards it being a technical error.

so "samples were contaminated," you say manchine that return the sequence do not make difference of DNA (Batteria) and RNA (HIV-1) virus??? possible error in data copy paste

HIV inserts its genome into the nuclear genome by using an enzyme, reverse transcriptase,  which copies the RNA as DNA. This same enzyme is used by humans to sequence RNA (the reverse transcribed version is called cDNA). Whereas Nanopore can actually read RNA the required processing for DNA samples is incompatible, so it was either gDNA or cDNA contamination during sample prep.

for me could be one lab escape... it seems the page is for me this https://science.sciencemag.org/content/266/5193/1981/tab-e-letters  the genbank file is LOCUS       WP_142766724             415 aa            linear   BCT 07-JAN-2020
DEFINITION  DDE-type integrase/transposase/recombinase, partial [Klebsiella
pneumoniae].

Bacteria are so numerous and so diverse that pretty much all sorts of oddities can happen, such as naked viral cDNA being absorbed and integrated via native recombinases. However, there are all sorts of "oh that's interesting" in environmental bacteria, but the truth is that only stuff that gets made into technology is looked at and that step requires a lot of effort. Say new enzymes, I have speculated the function of scores of enzymes and every time it I or others have placed then in the category "too much effort over publication impact".


text:
Consensus on how to select SNPs for creating a phylogeny when the full dataset is too large
I have got a quite large VCF file (about 20k samples with 160k SNPs after filtering for quality etc.) and I would like to get a phylogeny for it. However, the whole dataset's too large for my computational resources (unless I go for fasttree which would be my last resort). Therefore, I need to subset the SNPs and there's plenty of possible ways to do this (dump all rare SNPs with MAF below a threshold, only use SNPs in coding regions, only use missense SNPs; just to name a few).
Is there a consensus (or at least some comparative work or recent review; I found nothing) on what method works best and gives the tree that's closest to using the full dataset?

Yes, the SNPs have already been filtered on a handful of metrics.

OK - what species are you working on? It's worth looking at a paper that tries to build a tree using the same species as you and then following some of the steps they have performed. For example this paper if you are working with viruses.

I will give an answer that is specific to closely-related microbial genomes. Core-SNPs, SNPs which are within conserved and likely orthologous sequences, have been shown to be relatively informative for constructing phylogenies of large sequence sets[1].
Parsnp is the most popular tool for constructing core-genome alignments on microbial genomes, and extracts the core SNPs from the alignments. Parsnp then uses raxml or FastTree in order to produce a phylogeny from the core SNPs. I will say, though, that 20k sequences may be too many for Parsnp to run on (depending on how large the sequences are and your available resources).
Disclaimer: I am the maintainer of Parsnp (:

The subset that gives you the tree closest to the one you obtain with the full dataset does not necessarily is the tree nearest to reality... It could be that you get a more realistic result using only partial data for multiple reasons. Your best chance is in my opinion, as it usually is, making lots of tests and contrasting hypotheses.


text:
ATAC and RNA seq Up regulated region gene and scatter
I m not sure if I'm going the right way or not but when the gene and region and both upregulated in RNA seq and ATAC seq data I should be seeing  a positive correlation but it seems other wise.
Am i doing it in a wrong way or this can occur?

Adding the scatter with points

Updated Plot
But the issue is with the R^2 value it is very low

this looks like amazing results if true. How many data points do you have? can you add that to the plot or do a heatscatter ?

1400 UP region and genes ..I okay I will add scatter "this looks like amazing results if true." Im confused..

@StupidWolf updated the plot

this looks really weird.. did you flip the direction of the fold change.. and was the data subsetted somehow

Bit of background ..I did LRT test with HSC(Stem cell) as my reference level Then compared it to the rest of the Cell type which is Myeloid Progenitor,Leukemic stem cell, Monocyte ,Blast. So my idea was to see how the global changes happens across the normal and malignant cell types. So is the comparison is logical?

"did you flip the direction of the fold change." ..even if that happens would it change ?like if RNA fold change is changed to ATAC or vice versa?

if i go by this paper  figure 2 the fold changes are my figure looks similar except the slope  https://epigeneticsandchromatin.biomedcentral.com/articles/10.1186/s13072-019-0260-2/figures/2

For ATAC seq what i did is..if i have multiple regions annotated to gene i take the average of the region and then did the this fold change scatter

now looking at the dataframe looks like i made some blunder, " and was the data subsetted " this was the error

the last plot looks about right. most likely you need to subset on genes that are more highly expressed (because the range of the fold change looks crazy), and maybe try promoters for a start

"and maybe try promoters for a start" will do that


text:
How do I create a VCF file of all known pathogenic mutations in a gene of interest?
I would like to create a list (probably .vcf format would be good) of all known pathogenic missense mutations in a human gene of interest and then add other variants that could lead to the same pathogenic amino acid substitutions. I intend to search targeted sequencing data of my gene of interest for variants that would cause these particular substitions.
I'm very new to bioinformatics so I'm not so familiar with the appropriate databases and tools to use, but this is a rough idea of how it could be done:

extract all mutations listed in Ensembl within my gene of interest that are labelled with HGMD_MUTATION  (is there an option to download a vcf file with all HGMD mutations within my gene? I've only managed to do so for individual mutations.)
annotate list with reference and mutant amino acids for every mutation in the list
add additional mutations (which may not be present in HGMD) which could lead to the same substitution. E.g. I'm interested in the substitution Proline (reference codon: CCG) to Leucine (TTR or CTN) at AA position 102, so I would like to list all hypothetical mutations that would create a leucine codon in this position, not only the one given in Ensembl (which is CCG -> CTG).

Could you point me towards the best tools to use for these purposes? Unfortunately I haven't got a HGMD subscription. Many thanks!

Why HGMD specifically (not clear to me)? You can simply use the ClinVar database to find all pathogenic mutation within regions of interest/gene (you can check the genomic coordinates of any gene using either UCSC or ensemble).

I'm not sure how to generate the additional mutations, but I would say that HGMD is not the way to find all the pathogenic variants. I would probably filter this table by either Clinical significance or Evidence->Phenotype association.


text:
ambiguous calls in nanopolish
I am converting scores into frequencies (from nanopolish results) using Nanopore scripts.
The script filters out lots of rows on the premise that reads where LLR < T * N are ambiguous.
where:
LLR = log-likelihood ratio (log-likelihood methylated / log-likelihood unmethylated)
T   = arbitrary threshold. Above that value the C is methylated.
N   = number of Cs (CpG) in the fragment analyzed.

What? Why? Where?

What's the question here? If there is not enough support for a position, it gets filtered out...

I'm assuming that you are using the workflow defined here: https://nanopolish.readthedocs.io/en/latest/quickstart_call_methylation.html . Can you give more information about where you are in the process of the "Nanopore scripts?" I am guessing that you are at calculate_methylation_frequency.py based on available information, and in that case it is as simple as Wouter de Coster says.

Assuming that it is the postprocessing script calculate_methylation_frequency.py that you are using, user Wouter de Coster's comment is correct that the methylation event is filtered out at the frequency estimation stage if it has a low LLR. It will have  a low LLR if it is considered not confident enough to report. This means not only LLR > T * N, but also the statistical threshold LLR > 0 (see below).
Specifically, here is the code in the script:

def update_call_stats(key, num_called_cpg_sites, is_methylated, sequence):
    if key not in sites:
        sites[key] = SiteStats(num_called_cpg_sites, sequence)

    sites[key].num_reads += 1
    sites[key].called_sites += num_called_cpg_sites
    if is_methylated > 0:
        sites[key].called_sites_methylated += num_called_cpg_sites

That is the function that decides whether to add a site to output, which depends on is_methylated, which is set here:
    # [skip several lines of code]
        num_sites = int(record['num_motifs'])
        llr = float(record['log_lik_ratio'])

        # Skip ambiguous call
        if abs(llr) < args.call_threshold * num_sites:
            # [MP NOTE: SKIP if LLR < T * N]
            continue
        sequence = record['sequence']

        # [MP NOTE: SKIP IF LLR <= 0]
        is_methylated = llr > 0

        # [skip several lines of code]

        update_call_stats(key, num_sites, is_methylated, sequence)

It's a little messy to read because the sites dictionary is a global outside of update_call_stats().
So overall, it seems like if you want sites included that are below this threshold, you will have to either:

use an earlier stage of processing (output of nanopolish call-methylation I believe according to the docs) or
alter the script to change the heuristic thresholds.

Thanks! It wasn't easy to follow the script. I understood that the whole fragment is discarded if it can be defined which CpG is methylated (if at all). @Wouter de Coster, is that correct?


text:
TCGA dataset: different accession IDs mapped to same location?
I'm currently working on TCGA miRNA dataset. After constructing a reads matrix, I'm trying to find the isoform sequence. In my data, I have the genomic location (isoform_coords). I found that entries with different accession ID (MIMAT) were mapped to the same location.
I looked at the source files to make sure there were no errors during merging.
The data frame looks like this:
miRNA_ID        isoform_coords                      miRNA_region
hsa-mir-143     hg38:chr5:149428962-149428983:+     mature,MIMAT0000435
hsa-mir-143     hg38:chr5:149428962-149428983:+     mature,MIMAT0004599

I searched miRBase using the MIMAT ID and found that MIMAT0000435 and MIMAT0004599 are on different arms (5p and 3p) of hsa-mir-143:
hsa-miR-143-3p MIMAT0000435
hsa-miR-143-5p MIMAT0004599

This seems strange since there should be some distance between the 5p and 3p sequences.
Is there an explanation for this?
I need to get the sequence for each entry. I can do this using the Ensembl API. Is there a way I could perhaps do some mapping to get the 5p and 3p arms (so far I just did a manual search but I think there should be a more credible way to do this)? I would like to later analyse the expression patterns of the isoforms in cancer tissue so it would be useful to have isoform IDs as specific as possible.


text:
RNA_Seq Analysis in R, propmapped function issue
I'm currently trying to learn RNA-seq data analysis and differential expression in R. I've been using the Combine lab tutorial.
Everything was running well until the tutorial asked to use the propmapped function (on this page) on my BAM files and then view the promapped proportion for each sample. In this tutorial we should be getting promapped proportion ranging from 0.05 to 0.150, however in my case I have no reads that align to my mapped genome (promapped = 0 in my case, for each sample). I don't know why as I've been using the data they provided for this tutorial.
I was wondering if someone could follow the first few steps of this tutorial until the propmapped step and tell me if they eventually get the right proportion of reads or not after the propmapped function.
This could tell me if the tutorial is the culprit or not.

There must be a problem in the BAM files you are using. Here are steps I would do to figure where exactly is the problem.

Check the data you have downloaded. The tutorial is using 12 fastq.gz files (SRR1552444 - SRR1552455). You can check on sra how big they are supposed to be (For example SRR1552444 is supposed to be an ~1.7Gb file) and if your downloaded files are about the same. Check if the files are fine (In the terminal you can type zcat data/SRR1552444.fastq.gz | head to check the first few reads). If the files are 0 sized or somehow corrupted, you need to redownload them.

Check if your mapping reference looks alright (size wise, look at the first few lines) and if the reference is indexed (there is a .fai in the directory).

Check the size of the BAM files, they also need to be in the range of a few hundred Mb to a few Gb. I expect these will be corrupted for sure, otherwise, you would get your numbers. You might also want to check if they are indexed (although this should not be a problem, I suspect they will get indexed automatically). If there is nothing wrong the bamfiles, you might want to check another tool to summarise a bamfile (like samtools flagstat) to check if at the end the problem really is in propmapped. If that's the case.

This should at least help you to zoom to the core of your problem if not resolve it.

Thank you, i'll look into that


text:
Getting two class of promoters : ATAC seq data
In this paper they classify the promoters into two categories, highlighted in italics below.

Two Types of Promoters
Mammalian promoters can initiate transcription across broad or narrow
genomic regions, and these promoter shapes, broad or narrow (sharp),
correlate with distinct transcriptional regulatory mechanisms.33 We
evaluated cardiac promoters predicted from CAGE clusters for these 2
major types of promoters by calculating the interquartile range of
promoter CAGE clusters by determining the base pair distance between
10% and 90% of a promoter’s total signal. We observed the expected 2
distinct populations, defined as sharp (interquartile range <10 bp)
and broad promoters (interquartile range ≥10 bp; Figure 2A). Broad
promoters were those associated many different cellular functions,
including housekeeping functions. In contrast, genes with sharp
(narrow) promoters were those encoding proteins critical for
tissue-specific functions seen by the presence myofilament, muscle
tissue, and muscle contractions genes (Figure 2B). Thus,
tissue-specific genes important for LV specification and function were
more likely to have sharp promoters.

Now My question is can I do this only based on ATAC seq data. If yes how do i get that information of narrow and broad promoters.

Its highly unlikely that you will be able to categorize promoters on this basis using only ATAC data. This is a consequence of people using the terms "promoters" and "transcription start sites" interchangeably, when they are actaully different things. CAGE data measures transcription start sites which is the specific location where transcription starts. CAGE does not identify promoters, which are the complete region around the start of a gene that contains, at a minimum, all the sequences necessary to assemble the RNA polymerase II pre-initiation complex. Some definitions of promoter are even broader, including all regulatory sequence proximal the the gene start.
ATAC identifies regions of open chromatin, and I would expect the open chromatin regions around the starts of genes to be more like the broadest definition of promoter, and therefore to be significantly larger than the transcription starts sites they contain, whether they be broad or narrow TSSs.
You could of course look at a distribution of the peak sizes over TSSs and see if you found a bimodal distribution. If you did you could correlate the genes with broad/narrow TSS peaks with their known functions so see if it fitted the idea that narrows TSSs

thats a good start for me to start


text:
Reads count in metagenomics
Background:  I am developing a pipeline for metagenomic studies of human gut microbiote. In particular, I am mapping the reads data originating from shotgun whole genome sequencing onto a gene catalogue (similar to IGC) and counting reads mapped on each feature/gene. Till now I have been using NGLess, which is a rather convenient tool. However, we have some questions about how exactly it counts reads and it might be not the best option for scaling up our pipeline.
Question:
I am looking for alternative tools for gene count. A brief search has shown that htseq-count and featureCount are two frequently used tools. However they mainly appear in RNAseq context, and I am not certain whether they could be used in my case. I would also appreciate comments on the possibility of using samtools mpileup.

I think that there is a good argument for using technologies adapted from RNA-seq for metagenomic abundance quantification. The tool that I am most familiar with (that I use for this purpose) is kallisto, which uses pseudoalignment instead and is therefore less compute intensive.
For some background on the conceptual linkages between metagenomic abundance and transcriptomic abundance you can see the kallisto metagenome abundance paper or the explanatory blog post.
I think that it is preferable to use a tool that is specifically adapted for metagenomics, as there are always some numerical issues that might have been accounted for by the authors. I think that for this purpose it is maybe best not to naively use tools like mpileup, though I have no specific objection to it. However, if you look around you can see various papers using e.g. Cufflinks and Cuffdiff for effectively the same purpose, and Cufflinks is definitely intended for RNA-seq.
Depending on your specific application, you might want to use a different tool like mOTUS2, if you are interested in estimating abundance of specific genomes or lineages in a metagenome mixture. This method is probably more accurate as it is based on averaging across the whole genome, but it is not very useful if what you are interested in is more e.g. functional profiling of the metagenome.

Thank you for answering. We did try motus2, but it is restricted to a custom-made gene database and the marker genes grounded in the associeted species classification. What I am looking for is a tool for creating a gene count table, which we then could use to develop our own methods of species identification.

@Vadim In that case, I think that a tool like kallisto (or even cufflinks) might be most helpful. These will be providing counts as TPMs (though they also provide raw counts), of course, so if there is some statistical issue with that kind of compositional count then it might be a problem. However, metagenomics is intrinsically compositional so I'm not sure other methods will be much better. If you are ultimately going to do deconvolution of some kind, it might be helpful to look at what metabat2 is doing (more for ideas than for usefulness as a tool, as it's solving a different problem).