text:
Best practice for running GATK VQSR on X chromosome
According to GATK best practice, it is recommended that different VQSR models be built for SNPs and INDELs, because the annotations for high-quality SNPs and INDELs are systematically different (if I understand it correctly). Since annotations for good variants on autosomes could be different from those on X chromosome, e.g., DP for good variants on X chromosome could be substantially smaller than DP on autosomes due to having large number of male samples, it seems reasonable to build VQSR models separately for X chromosome and autosomes. However, no such advice is proposed on GATK website.
My question is:
Should we build VQSR models separately for autosomes and X chromosome? If yes, we will have 4 VQSR models: Auto SNPs, Auto INDELs, X SNPs and X INDELs.
Cross-posting: https://www.biostars.org/p/9483110/

I see no reason to treat different as long as all samples are derived from the same sex in  humans or have the same ploidy generally speaking. I agree males differ in DP values, which are about halved compared to females, and this might impact other covariates, so for X chromosome I would do VQSR on only males, and in separate for only females. But I would not make any distinction between autosomes and X for females,

@JRodrigoF Could you elaborate on your workflow? Since VQSR is a site filtering procedure, once you have joint-called X chromosome with both male and female samples, there is no meaningful way to "separate" them for VQSR.

From the very beginning of the GATK workflow you have to make a difference between females and males if you ought to be inclusive of the X chromosome. Males have to go through haplotycaller using the argument ploidy = 1 , whereas females can use default value = 2 , ... then from there you already have them in separate , even if you joint call them altogether at the genotypeGVCF level, you always know which of your samples are females and which males, right? So of course you can separate them again into two different multi-sample VCFs using bcftools for instance,

@JRodrigoF Thanks for your comment. What I meant in the last comment is that, if I joint call male and female samples together, there will be only one value per site for each site annotation that is considered by VQSR. Even if I separate male from female samples, the site annotations will not change except trivial ones such as AC, AF, etc.. Did you suggest joint calling X chromosomes for males and females separately and combining them after VQSR? If you feel confident about your comments, would you consider converting it to an answer?

Yes. Essentially boils down to treat X chromosomes from males in separate and from the very beginning, that is at the haplotypecaller step, then you follow by joint calling (genotypeGVCF), then VQSR (SNPs and INDELs) , and only after this merge with the rest of the autosomes from both males and females and with the X chromosomes from females,
OK. Will turn it into an answer,

Have you had a look at the thousand genomes papers? The VQSR was used for that dataset, including on the X chromosome (PAR1, PAR2, and all). The supplementary methods sections are extensive.

@winni2k Did you mean the 2015 1000G phase 3 paper? Or newer technical updates to 1000G? (GRCh38 realigned version, high-coverage version, etc.)

I mean the phase 1 and phase 3 papers from way back ;) VQSR is old. This should not be a new problem :)

You are correct that differences in ploidy between autosomes (always ploidy=2) and X chromosome mean that the same VQSR model cannot be applied to both cases equally.
However, this is the case only for male X chromosomes (in human species).
Therefore, yes, I'd recommend building a separate VQSR model when it comes to male X chromosomes. Importantly, the processing has to be different not only at the VQSR step, but from the very beginning at the HaplotypeCaller step (if you are using GATK also for variant calling), which has to be used with the option ploidy = 1.
So,
On one hand -> Processing of autosomes (females and males) and x chromosomes (females only) can be done altogether with no differences all the way to VQSR (SNPs + INDELs as in best practices)
On the other hand -> Processing of X chromosomes (only males) goes in separate, from variant calling to VQSR (SNPs + INDELs as in best practices)

If you're going to go for this, do remember that you will need to break up the X chromosome into three sections, as there are two pseudoautosomal regions PAR1 and PAR2 at the beginning and the end. See https://en.wikipedia.org/wiki/Pseudoautosomal_region


text:
GWAS : input files for rrBLUP package in R
I want to do a verification of the result of an RNAseq analysis of human cancer samples by GWAS, I want to use rrBLUP package in R, so I need three input file genotype, phenotype and mapping file, so I am confused about how I can get these files.

source: Stack Exchange

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

sort: votes


text:
Keeping rows with a cut-off in all columns
I have a data frame of raw read counts, genes in rows and samples in columns
> dput(head(oc))
structure(list(sample1 = c(13L, 0L, 1L, 42L, 15L, 1137L), sample2 = c(21L,
2L, 20L, 51L, 5L, 203L), sample3 = c(25L, 2L, 6L, 71L, 36L, 1346L
), sample4 = c(59L, 3L, 8L, 58L, 42L, 3428L), sample5 = c(19L,
4L, 8L, 52L, 1L, 201L)), row.names = c("WASH7P", "MIR6859-1",
"CICP27", "WASH9P", "MTND1P23", "MTND2P28"), class = "data.frame")
>

> head(oc)
          sample1 sample2 sample3 sample4 sample5
WASH7P         13      21      25      59      19
MIR6859-1       0       2       2       3       4
CICP27          1      20       6       8       8
WASH9P         42      51      71      58      52
MTND1P23       15       5      36      42       1
MTND2P28     1137     203    1346    3428     201
>

I just want to keep genes with at least 10 reads in each sample, I mean even if a gene has less than 10 reads in one sample and the rest of four samples have 10 reads, this gene should be removed. For instance here the second, third and fifth rows should be removed
I tied this but all I obtained zero
> df=t(oc)
> cutoff <- 10
>
> # encode the data based on the cutoff
> df[df >= cutoff] <- 1
> df[df < cutoff] <- 0
> View(df)
> rowSums(df)
sample1 sample2 sample3 sample4 sample5
      0       0       0       0       0
>

> apply(df >= 5, 1, sum)
sample1 sample2 sample3 sample4 sample5
      0       0       0       0       0
> apply(df >= 1, 1, sum)
sample1 sample2 sample3 sample4 sample5
      0       0       0       0       0
>

Any suggestion please?

apply(my_df >= 10, 1, sum) would give you that how many of the samples satisfy your criteria. Basically you apply() the function sum() to the rows (denoted as 1). The sum() does the trick that @Devon Ryan suggested above.

I get all zero this is weird  > apply(df >= 5, 1, sum)
sample1 sample2 sample3 sample4 sample5
0       0       0       0       0

apply(df >= 1, 1, sum)
sample1 sample2 sample3 sample4 sample5
0       0       0       0       0

Check your df, looks like you have manipulated that data frame.

No just df=t(oc) I used in the code

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


text:
Find corresponding symbol for gene used in Seurat
I have a Seurat object
> pbmc
An object of class Seurat
48608 features across 11709 samples within 2 assays
Active assay: SCT (22220 features, 3000 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap
>

I want to plot TCRVδ2
> FeaturePlot(pbmc,"TCRVδ2")
Error: None of the requested features were found: TCRVδ2 in slot data
In addition: Warning message:
In FetchData(object = object, vars = c(dims, "ident", features),  :
  The following requested variables were not found: TCRVδ2
>

This data base is about human TCR genes
https://www.genenames.org/data/genegroup/#!/group/370
But I even can not find this gene here
Do you know what is the corresponding symbol for this gene used in Seurat?

TCRVδ2 is not a gene, but it defines a specific subset of gamma-delta TCRs. You can look at TRDC (delta chain) and TRGC1 and TRGC2 (gamma chains). TCRVδ2 should have the combo TRDC-TRGC1, while TCRVδ1 should express TRDC-TRGC2. You may need to increase the resolution of your clustering, to avoid clustering the gamma-delta T cells with the NK cells.

Thank you so much. I have another question about gene symbols. CD45 is corresponding to PTPRC butler most of CD genes I can not find the corresponding official gene symbols. Please, do you know how I can find corresponding gene symbols for CD genes?


text:
Avoiding the warning 'The following features were omitted as they were not found in the scale.data slot for the SCT assay' in Seurat
I have a Seurat object
I want to plot a list of genes but I got a warning saying most of genes removed
I have tried with both my genes of interest and all variable genes but the results is very similar
Does this mean these genes are basically are not expressed in my data?
   > t$V1
  [1] "IGH"          "IGHV1-2"      "IGHV1-3"      "IGHV1-12"     "IGHV1-8"
  [6] "IGHV1-14"     "IGHV1-17"     "IGHV1-18"     "IGHV1-24"     "IGHV1-45"
 [11] "IGHV1-46"     "IGHV1-58"     "IGHV1-67"     "IGHV1-68"     "IGHV1-69"
 [16] "IGHV1-69D"    "IGHV1-38-4"   "IGHV1-69-2"   "IGHV2-5"      "IGHV2-10"
 [21] "IGHV2-26"     "IGHV2-70"     "IGHV2-70D"    "IGHV3-6"      "IGHV3-7"
 [26] "IGHV3-9"      "IGHV3-11"     "IGHV3-13"     "IGHV3-15"     "IGHV3-16"
 [31] "IGHV3-19"     "IGHV3-20"     "IGHV3-21"     "IGHV3-22"     "IGHV3-23"
 [36] "IGHV3-25"     "IGHV3-29"     "IGHV3-30"     "IGHV3-30-2"   "IGHV3-30-3"
 [41] "IGHV3-30-5"   "IGHV3-32"     "IGHV3-33"     "IGHV3-33-2"   "IGHV3-35"
 [46] "IGHV3-36"     "IGHV3-37"     "IGHV3-38"     "IGHV3-41"     "IGHV3-42"
 [51] "IGHV3-43"     "IGHV3-43D"    "IGHV3-47"     "IGHV3-48"     "IGHV3-49"
 [56] "IGHV3-50"     "IGHV3-52"     "IGHV3-53"     "IGHV3-54"     "IGHV3-57"
 [61] "IGHV3-60"     "IGHV3-62"     "IGHV3-63"     "IGHV3-64"     "IGHV3-65"
 [66] "IGHV3-66"     "IGHV3-71"     "IGHV3-72"     "IGHV3-73"     "IGHV3-74"
 [71] "IGHV3-75"     "IGHV3-76"     "IGHV3-79"     "IGHV3-38-3"   "IGHV3-69-1"
 [76] "IGHV4-4"      "IGHV4-28"     "IGHV4-30-1"   "IGHV4-30-2"   "IGHV4-30-4"
 [81] "IGHV4-31"     "IGHV4-34"     "IGHV4-39"     "IGHV4-55"     "IGHV4-59"
 [86] "IGHV4-61"     "IGHV4-80"     "IGHV4-38-2"   "IGHV5-51"     "IGHV5-78"
 [91] "IGHV5-10-1"   "IGHV6-1"      "IGHV7-4-1"    "IGHV7-40"     "IGHV7-27"
 [96] "IGHV7-56"     "IGHV7-34-1"   "IGHV7-81"     "IGHVII-1-1"   "IGHVII-15-1"
[101] "IGHVII-20-1"  "IGHVII-22-1"  "IGHVII-26-2"  "IGHVII-28-1"  "IGHVII-30-1"
[106] "IGHVII-30-21" "IGHVII-31-1"  "IGHVII-33-1"  "IGHVII-40-1"  "IGHVII-43-1"
[111] "IGHVII-44-2"  "IGHVII-46-1"  "IGHVII-49-1"  "IGHVII-51-2"  "IGHVII-53-1"
[116] "IGHVII-60-1"  "IGHVII-62-1"  "IGHVII-65-1"  "IGHVII-67-1"  "IGHVII-74-1"
[121] "IGHVII-78-1"  "IGHVIII-2-1"  "IGHVIII-5-1"  "IGHVIII-5-2"  "IGHVIII-11-1"
[126] "IGHVIII-13-1" "IGHVIII-16-1" "IGHVIII-22-2" "IGHVIII-25-1" "IGHVIII-26-1"
[131] "IGHVIII-38-1" "IGHVIII-44"   "IGHVIII-47-1" "IGHV8-51-1"   "IGHVIII-67-2"
[136] "IGHVIII-67-3" "IGHVIII-67-4" "IGHVIII-76-1" "IGHVIII-82"   "IGHVIV-44-1"
[141] "IGHD1-1"      "IGHD1-7"      "IGHD1-14"     "IGHD1-20"     "IGHD1-26"
[146] "IGHD2-2"      "IGHD2-8"      "IGHD2-15"     "IGHD2-21"     "IGHD3-3"
[151] "IGHD3-9"      "IGHD3-10"     "IGHD3-16"     "IGHD3-22"     "IGHD4-4"
[156] "IGHD4-11"     "IGHD4-17"     "IGHD4-23"     "IGHD5-5"      "IGHD5-12"
[161] "IGHD5-18"     "IGHD5-24"     "IGHD6-6"      "IGHD6-13"     "IGHD6-19"
[166] "IGHD6-25"     "IGHD7-27"     "IGHJ1"        "IGHJ1P"       "IGHJ2"
[171] "IGHJ2P"       "IGHJ3"        "IGHJ3P"       "IGHJ4"        "IGHJ5"
[176] "IGHJ6"        "IGHA1"        "IGHA2"        "IGHG1"        "IGHG2"
[181] "IGHG3"        "IGHG4"        "IGHGP"        "IGHD"         "IGHE"
[186] "IGHEP1"       "IGHM"
> pbmc=ScaleData(
+     pbmc,
+     features = t$V1)
Centering and scaling data matrix
  |================================================================================| 100%
> DoHeatmap(pbmc,t$V1,group.by = "Condition")
Warning message:
In DoHeatmap(pbmc, t$V1, group.by = "Condition") :
  The following features were omitted as they were not found in the scale.data slot for the SCT assay: IGHEP1, IGHE, IGHJ6, IGHJ5, IGHJ4, IGHJ3P, IGHJ3, IGHJ2P, IGHJ2, IGHJ1P, IGHJ1, IGHD7-27, IGHD6-25, IGHD6-19, IGHD6-13, IGHD6-6, IGHD5-24, IGHD5-18, IGHD5-12, IGHD5-5, IGHD4-23, IGHD4-17, IGHD4-11, IGHD4-4, IGHD3-22, IGHD3-16, IGHD3-10, IGHD3-9, IGHD3-3, IGHD2-21, IGHD2-15, IGHD2-8, IGHD2-2, IGHD1-26, IGHD1-20, IGHD1-14, IGHD1-7, IGHD1-1, IGHVIV-44-1, IGHVIII-82, IGHVIII-76-1, IGHVIII-67-4, IGHVIII-67-3, IGHVIII-67-2, IGHV8-51-1, IGHVIII-47-1, IGHVIII-44, IGHVIII-38-1, IGHVIII-26-1, IGHVIII-25-1, IGHVIII-22-2, IGHVIII-16-1, IGHVIII-13-1, IGHVIII-11-1, IGHVIII-5-2, IGHVIII-5-1, IGHVIII-2-1, IGHVII-78-1, IGHVII-74-1, IGHVII-67-1, IGHVII-65-1, IGHVII-62-1, IGHVII-60-1, IGHVII-53-1, IGHVII-51-2, IGHVII-49-1, IGHVII-46-1, IGHVII-44-2, IGHVII-43-1, IGHVII-40-1, IGHVII-33-1, IGHVII-31-1, IGHVII-30-21, IGHVII-30-1, IGHVII-28-1, IGHVII-26-2, IGHVII-22-1, IGHVII-20-1, IGHVII-15-1, IGHVII-1-1, IGH [... truncated]
> pbmc=ScaleData(
+     pbmc,
+     features =rownames(pbmc))
Centering and scaling data matrix
  |================================================================================| 100%
> rm(pbmc)
> pbmc=readRDS("pbmc.rds")
> pbmc=ScaleData(
+     pbmc,
+     features =rownames(pbmc))
Centering and scaling data matrix
  |================================================================================| 100%
> DoHeatmap(pbmc,t$V1,group.by = "Condition")
Warning message:
In DoHeatmap(pbmc, t$V1, group.by = "Condition") :
  The following features were omitted as they were not found in the scale.data slot for the SCT assay: IGHEP1, IGHE, IGHJ6, IGHJ5, IGHJ4, IGHJ3P, IGHJ3, IGHJ2P, IGHJ2, IGHJ1P, IGHJ1, IGHD7-27, IGHD6-25, IGHD6-19, IGHD6-13, IGHD6-6, IGHD5-24, IGHD5-18, IGHD5-12, IGHD5-5, IGHD4-23, IGHD4-17, IGHD4-11, IGHD4-4, IGHD3-22, IGHD3-16, IGHD3-10, IGHD3-9, IGHD3-3, IGHD2-21, IGHD2-15, IGHD2-8, IGHD2-2, IGHD1-26, IGHD1-20, IGHD1-14, IGHD1-7, IGHD1-1, IGHVIV-44-1, IGHVIII-82, IGHVIII-76-1, IGHVIII-67-4, IGHVIII-67-3, IGHVIII-67-2, IGHV8-51-1, IGHVIII-47-1, IGHVIII-44, IGHVIII-38-1, IGHVIII-26-1, IGHVIII-25-1, IGHVIII-22-2, IGHVIII-16-1, IGHVIII-13-1, IGHVIII-11-1, IGHVIII-5-2, IGHVIII-5-1, IGHVIII-2-1, IGHVII-78-1, IGHVII-74-1, IGHVII-67-1, IGHVII-65-1, IGHVII-62-1, IGHVII-60-1, IGHVII-53-1, IGHVII-51-2, IGHVII-49-1, IGHVII-46-1, IGHVII-44-2, IGHVII-43-1, IGHVII-40-1, IGHVII-33-1, IGHVII-31-1, IGHVII-30-21, IGHVII-30-1, IGHVII-28-1, IGHVII-26-2, IGHVII-22-1, IGHVII-20-1, IGHVII-15-1, IGHVII-1-1, IGH [... truncated]
>

Specify your genes of interest using the features parameter of the ScaleData() function. If this parameter is not specified, only "variable genes" are scaled and that makes sense, genes that do not demonstrate variability across conditions is not interesting in general. Anyhow, since the  DoHeatmap() would be looking at the "scaled data slot", it will not be able to find your genes of interest unless they show variation.

Thank you. I use sctranscform, I should scale my data before or after sctransform function?

From their vignette: Note that this single command replaces NormalizeData(), ScaleData(), and FindVariableFeatures().

I checked sctransform with or without scale data and the results was different

1- Please do not edit your questions extensively, these questions meant to help future users too, now the answer and the question does not match. 2- Please read documentation associated with the functions/packages you use: SCTransform() doc says return.only.var.genes:  If set to TRUE the scale.data matrices in output assay are subset to contain only the variable genes; default is TRUE. 3- Please show some effort, you will learn and develop faster if you learn by doing instead of asking questions on each step.


text:
Python script to simultaneously generate multiple pdbqt files for AutoDockTools?
I have a folder of pdb ligand structures that I'd like to test in some docking experiments with a certain protein using AutoDock Vina.
I am not familiar with Python so I've always been using the interactive procedure on AutoDockTools to prepare the pdbqt ligand files instead of the prepare_ligand4.py script.
Could you please help me understand how to iterate this python script in a way that it can be applied to a set of pdb ligand files in order to create the corresponding set of pdbqt ligand files?
Thank you!

use Raccoon or Raccoon2 : https://autodock.scripps.edu/resources/raccoon/ , https://autodock.scripps.edu/resources/raccoon2/

Answer from @matteo-ferla, converted from comment:
If command line interface is more suitable for you, you may consider OpenBabel, which can predict Gasteiger charges, protonate at pH 7, accept a smi file of multiple SMILES or a multi-entry sdf file and can output in PDQBT format, but is not quite the same as the python script allows merging and various settings.
But to use the script iteratively (a horrendous thing to do as it loads up the python kernel each time) check out how to loop files via bash


text:
Creating new seurat object with new matrix
I'm trying to do a cross-species comparison between the patient TME vs. the PTX TME to understand gene expression conservation. I have already converted between mouse and human genes giving me a resultant gene expression matrix that is just one set of (homologous) genes X cells for all mouse and human cells.
musGenes <- rownames(mouse_only)
# Basic function to convert mouse to human gene names
convertMouseGeneList <- function(x){
require("biomaRt")
h = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
m = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
genesV2 = getLDS(attributes = c("mgi_symbol"), filters = "mgi_symbol", values = x , mart = m, attributesL = c("hgnc_symbol"), martL = h, uniqueRows=T)
humanx <- unique(genesV2[, 2])
# Print the first 6 genes found to the screen
print(head(humanx))
return(humanx)
}
genes <- convertMouseGeneList(musGenes)

Now how do I create a new Seurat object with the converted mouse gene list, such as rename the matrix, so that the rownames/colnames are human symbols, rather than mouse symbols?

is this thread useful? They do rename rownames of the just created seurat object.


text:
R limma alternatives in Python
The R package limma is ideal to perform differential expression analysis. Is there any limma alternative in Python?
I'm trying to use statsmodels and scikitlearn in conjunction with some other cool tools (such as pycombat) to get limma-like workflows. However, I'm missing many capabilities from the limma toolbox.

Is there a possibility for you to call R from Python? If so you can try Rpy2 (https://pypi.org/project/rpy2/)

That's true. You can always wrap Python around R. For example, if working with limma, you can use Python to prepare the pieces that limma needs to create an Expression Set object with the ExpressionSet() function (expression, phenotype and feature data). However you will be always forced to write the limma workflow using R syntax. In any case I will give it a try! Many Thanks!

I am working on EXACTLY this thing. My methylize python package supports some of these functions, but I found this while looking for packages that might give me idea about what features people want.

August 2021
To date, it seems that the response to this question is "No". However, there is a GitHub repository called edgePy, aiming to "become an implementation of edgeR for differential expression analysis in the Python language" (edgeR is somewhat similar to limma).
PS: Many thanks to @ATpoint90 (the Twitter "hive mind" summoner) and @sridnona6 (edgePy developer).


text:
sctransform - mitochondrial expression filtering / transformation
I have a Seurat object which has a high expression of mitochondrial genes
> pbmc1
An object of class Seurat
36601 features across 18338 samples within 1 assay
Active assay: RNA (36601 features, 0 variable features)
>

I want to use SCTransform in Seurat, I don't know if still I should define the mitochondrial percentage or not, if not, how SCTransform function knows which percentage suits for my data?
I used SCTransform without clipping MT genes and I noticed I have the same number of cells
> pbmc_SCTransform <- SCTransform(pbmc, method = "glmGamPoi", vars.to.regress = "percent.mt", verbose = FALSE)
> pbmc_SCTransform
An object of class Seurat
58629 features across 18338 samples within 2 assays
Active assay: SCT (22028 features, 3000 variable features)
 1 other assay present: RNA

Is it fine or I should clip MT genes before running SCTransform function?
I am not too sure about the input composition but from 10X websummary two of my samples likely suffer from poor sample prep or poor beginning cell health. The gene expression median genes per cell are 41 and 50 and differential expression show only mitochondrial (MT) and MALAT1 genes. This is indicative of subpar lysis conditions.

Frustrating, why down voting? Do you know the answer? You believe in Seurat vignette  this has been already mentioned? No

Mitochondrial gene filtering is done on a case-by-case basis, which is why there's a bit of discussion about choosing appropriate thresholds in the Seurat PBMC3k vignette.
In your specific case, there are a lot of cells with very high mitochondrial gene content. It would be a good idea to work out why that is before doing any more filtering.
However, you're asking specifically about sctransform, and the SCtransform vignette does account for mitochondrial gene content in the transformation:

During normalization, we can also remove confounding sources of variation, for example, mitochondrial mapping percentage

# store mitochondrial percentage in object meta data
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")

# run sctransform
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)

There's no filtering in that vignette, but there's also no mention of how clean the data was initially. It's still a good idea to look at your data for outliers, try to understand why they are outliers, and have a think about whether or not to exclude them.


text:
Barplot without alphabetizing x-axis
I'm inputting a .csv file and then plotting a barplot based on the colnames but when I go plot the order the names are being alphabetized rather than get the order that is in the .csv file. How to I prevent that from happening?
p2 = ggplot(data=clusters_per_sample, aes(x=sample, y=Perc, fill=cluster)) +
  geom_bar(stat="identity")+
  scale_fill_manual(values = getPalette(22))+
  theme_minimal() + labs(y="percentage of total # cells in the sample") + theme(axis.text.x = element_text(angle = 90, face = "bold")) + theme(axis.text.x = element_text(size=12, color="black"))

I tried this to get rid of the alphabetizing but my bars are very thin when doing so:
p2 = ggplot(data=clusters_per_sample_2, aes(x=sample, y=Perc, fill=cluster)) +
geom_bar(stat="identity")+ scale_x_discrete(limits=clusters_per_sample_2$sample) +
heme_minimal() + labs(y="percentage of total # cells in the sample") + theme(axis.text.x = element_text(angle = 90, face = "bold")) + theme(axis.text.x = element_text(size=12, color="black"))

Any suggestions?

ggplot2 uses factors on the underlying data frame. To have the data plotted in a custom order, apply a custom ordered factor on the underlying dataset.
For example, if clusters_per_sample$sample is your custom-ordered vector, do this:
clusters_per_sample$sample <- factor(clusters_per_sample$sample, levels = unique(clusters_per_sample$sample), ordered = TRUE)

Now when you plot sample as an aesthetic, ggplot has a factor component to work with already, so it won't try and make a new one (this new factor would have been in alphabetical order, which we circumvent with the ordered=TRUE parameter in our custom factor).

Where do I add that line in my code?

Add it in before passing the object to ggplot(). This step is part of preparing your data for plotting.


text:
Hisat2 compatibility for long reads (Pacbio)
I am working with a high performance cluster computer containing 112 threads. I am trying to align PacBio transcriptome reads against the genome to count the gene number. For pair end read i used the following workflow:
#convert gff to gtf
/home/software/cufflinks-2.2.1/gffread xxx.gff -T -o xxx.gtf
#build index
/home/software/hisat2-2.2.1/hisat2_extract_exons.py xxx.gtf > xxx.exon
/home/software/hisat2-2.2.1/hisat2_extract_splice_sites.py xxx.gtf > xxx.ss
/home/software/hisat2-2.2.1/hisat2-build -p 20 xxx.fa --ss xxx.ss --exon xxx.exon xxx_tran
#align reads
/home/software/hisat2-2.2.1/hisat2 -p 80 -x ./xxx_tran -1 R1.fq -2 R2.fq -S gonad.sam

#convert to bam and sort by read name
samtools view -bS gonad.sam -o gonad.bam

samtools sort -@ 80 gonad.bam -o gonad.sorted.bam

#count gene
htseq-count -f bam --strand=no gonad.sorted.bam xxx.all.gtf > gonad-counts.txt

But in case of Long read i am not sure what i need to do during the alignment stage. I have tried with minimap2 like as following:
#align reads
/home/software/minimap2-2.17_x64-linux/minimap2 -ax splice:hq -uf yyy.fa gill.fq > gill.sam

in following, i stucked in counting stage:
#count gene
htseq-count -f bam --strand=no gill.sorted.bam yyy.gtf > gill-counts.txt
error txt:
[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
438395 GFF lines processed.
[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'
Error occured when processing input (record #79 in file gill.sorted.bam):
  'NoneType' object has no attribute 'encode'
  [Exception type: AttributeError, raised in _HTSeq.pyx:1379]

My question is that during counting it requires sorted bam and genome gtf file that was generated succesfully during previous staps, why it is getting failed even though the same pipeline is absolutely fine for pair end reads? Hisat2 can align long read? Is there any other way do that? As a novice i will be always grateful to you. Thanks.
Update: I tried generating an index:
# generate index
samtools index gill.sorted.bam

And the output was a .bai file like "gill.sorted.bam.bai". And I tried again to count as written above. But then I found the following error:
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
438395 GFF lines processed.
Error occured when processing input (record #79 in file gill.sorted.bam):
  'NoneType' object has no attribute 'encode'
  [Exception type: AttributeError, raised in _HTSeq.pyx:1379]

Is there any way?
Update 2:
Here is what the BAM file looks like:
SRR6418045.89565    2048    PYXC01000001.1  9893    60  5794H11M1I42M1I3M2I19M1D20M1I15M1D20M3D43M1I11M3D8M1I19M3I24M1I7M1I4M1I5M1D6M1D14M1I45M1D2M1D3M1D4M1I8M1I25M1D21M1D62M1D3M1D25M2D21M2I30M1I25M1D2M1D5M2I4M1D13M1D12M1D10M1I59M21564H    *   0   0   GAATACATAAACTGGATAGCGCACTGTTAATATTTACAGTGAAGAGAGTCTCTTCGAATTTTCCTCCTCAATCATCTTATTAGAAAGATGTTGAGTGGGCAGCAGCCAAATAAAATACTGTTCAGTGATCACTTCCAATGATGCTGTGCACTTACAGAAGCATGTATGGTACACAGTGTTATGCCTAATTTAAAAAACTAATGTAAATTATTTAAAAGTTTTTTTTTTTTTTATGAACTTTTTAGCCACATTGAATGGCTTCAGAATATCTACAGCCTAATGTATTTTCTTAGTAGCCTAATGAAGAAAGAATATTACCTTAAGTTACCATAGATTGACGCCAAAACTCGTCAAATCTCAAACCCCTAACATTTGTTCAACTGTTGTATACATGCTCATACACAGTGTTTCTTGTGTTCCTTAACCCACCGGGTAAAAAGTGGAGGTCCATGTCTGATACGTAGTTATCAATGTACTGTTGCCCATTTAGAGAGCACTTTGGTTTTCCCTCACCAGATAAGGTTTTTTCTCACTTTCAGTGTAAACCTGAGACTTGTTTGCACTTAATATTGGGGATTGTGTAACTCTTATTCAAGGAATTGCCTGTTTTCCAATTTGGTCTGATTACACTAACCTTTGTGTAAAATTATGGTTCAAATAAACTGTATCACA    -)/////,+$/"/..+//-/-+/.-/-.*,./.,)//'//.+/./..-/./../"/*/.),#)*.)-+)*-(///&.//.*'/,*////*-,./+.++#",.///,/,*,/-,)/,/)/.-)/.///./-..-*++)-(/*,../..)+//.////+*...////./...(././,/.,*.-'*/)(.-*)&(#'&,(#'....').,,$*')+#&,,)%%%$&,)*)*(*,..,'/-.+*)*,)+',-).-)*//,-",-+./)///'-.//,-+,.*//.-/%,**)-..../.*/,*'/.*//,*/.,(&.,+-)-'-*/*-/.-/-&+$&,'.'.&,.*,)/-."/,.,*/-//.,)/*,%.+*)/+.,--,./*--/.-.//'///)/--/$&//-//.//.-)/),,./-)+&../*,*,..-/+,)/,,,-*,/-,/.(.(///./(////./&-./-/+//*...////,&%+$&%#%$$--.././/.-+',-,#(.&%(-/,*/.///*..--++*)/-.*..*)*(,$*/--*.+./'//...-..)-.&.-)&))+.-,&#'-'--..'/'-//',//-*.)/+-*....%../.'(/*-$./,(.(////.#/..+/+%/'-'%...//*$+#--+/&-,-/,*-.-*/-.&../-,%/    NM:i:55 ms:i:343    AS:i:343    nn:i:0  ts:A:+  tp:A:P  cm:i:77 s1:i:460    s2:i:0  de:f:0.0656 SA:Z:PYXC01000001.1,9893,+,4147S672M2I23209S,60,37;PYXC01000001.1,9893,+,847S669M26I26488S,60,42;PYXC01000001.1,9893,+,26975S673M3I379S,60,44;PYXC01000001.1,9910,-,10980S654M9I16387S,60,41;PYXC01000001.1,9900,-,20736S661M9I6624S,60,43;PYXC01000001.1,9893,+,7409S673M15I19933S,60,49;PYXC01000001.1,9893,-,27284S667M12I67S,60,52;PYXC01000001.1,9893,-,19099S673M3I8255S,60,48;PYXC01000001.1,9893,-,7656S673M1652I18049S,4,1688;PYXC01000001.1,9908,+,2517S649M17I24847S,60,51;PYXC01000001.1,9893,+,23715S672M1D3643S,60,51;PYXC01000001.1,9893,-,25633S669M22I1706S,60,53;PYXC01000001.1,9893,-,17479S668M1I9882S,60,53;PYXC01000001.1,9893,-,14214S673M25I13118S,22,55;PYXC01000001.1,9899,+,13933S663M13I13421S,60,55;PYXC01000001.1,9900,-,22352S658M21I4999S,60,55;PYXC01000001.1,9901,+,15556S665M7I11802S,60,60;PYXC01000001.1,9893,+,12284S668M9I15069S,60,58;PYXC01000001.1,9909,-,24008S657M16I3349S,60,57;PYXC01000001.1,9902,+,25358S661M10I2001S,60,60;PYXC01000001.1,9893,+,10659S671M2D16700S,60,60;PYXC01000001.1,9893,+,9039S668M10I18313S,60,51;PYXC01000001.1,9893,-,2790S673M21I24546S,60,65;PYXC01000001.1,9899,+,22114S669M8I5239S,60,59;PYXC01000001.1,9893,-,6025S673M17I21315S,60,64;PYXC01000001.1,9908,-,12604S645M13I14768S,60,64;PYXC01000001.1,9893,-,4418S673M3I22936S,60,71;PYXC01000001.1,9893,-,1162S673M15I26180S,60,70;PYXC01000001.1,9898,+,17178S641M1I10210S,60,63;PYXC01000001.1,9893,+,20477S673M13I6867S,60,67;PYXC01000001.1,10390,-,64S160M2D27806S,60,8;   rl:i:2011
SRR6418045.89565    2048    PYXC01000001.1  9893    60  12284H12M1D27M1D9M1I15M1I26M1D19M2I16M2I2M1I17M1I32M1I31M2D8M3I25M1I6M1I19M2I12M1D32M2D15M1I6M2I7M1D4M2D17M1I10M1I32M1I36M1D4M1D2M1D20M1I59M1I11M1D10M1D10M1I24M1I9M2D2M1D28M1I24M1I11M15069H   *   0   0   GAATACATAAATGATAGCGCACTGTTAATATTTACAGTGAGAGAGTCTACTTGAATTCCTCCTCAAATCATCTTATTTAGAAAGATGTTGATGGGAGCAACCAAATAAAAATCTACTGTTCAGTGATCCAACTTTCACAAAATGATGCTGTTGCACTTACAGAAGCATGTATGGTACACAGTGTTATGACTATTTTAATAAAAAACAATGTAAATTTTAAAAGGTTTTTTTTTTTTTATGAACTTTTTAGCAACATTGAATGGTTCAGAAATATCTTAGCAAGCCTAATGTATTCTTAGTAGCCTAATGAAGAAAGAATATTACTAAGTTACCAGTAGGAATCTGAACGCCAAAATGTCATCTCAAACCCCTAACATTTTGGTTCAATCTGTTGTATACATGCATCATACACAGTGTTTCTTTGTGTTCCTTTAACCACCGGGTAAAAAGTGGAGGTCATGCTATTACTGTAGTTATCAATGTAACTGTTGCCCACTTTTAGAGAGCACTTTGGTTTTCTCACCAGATAAGGTTTTTTCTCACTTTTTATGTAAACTGAGACTTGTTGCACTGTACTATATTGAATATGTGTAACTCTTATTTTCAAGGATGCCTGTTTTCCATTTGGTCTATTAACACTAACACTTTGTGTAAACTTATGGTTCAAAATAAACTGT   -*//..//-*./%/,/./,,///.,/+//*/,////..,*..+/..-,(*.,.%&".*./%-..((//$///...)+-&.+'/////.+*$(),,../$.'.,)//-++$.."//,-.,/////'//,/$...*//./-()./+//,///,...%/-//,*.*/.//&/.//...///./..(-.///#/.+$(#.'..((%*(.)(-+-..%.*,((,,(..(-&+*)*******..+/$-..,++..--$-..'-.)//..-)/.,,*...//./.#.&/)*/)*(/./*,.&-(+////*/.).).+//-*//*)/*.*..-(//-/&,/-//..)*//.+.$/.',-+,(+(&//..//.,).%')./*///.--)--*.'(/"//...//.././/-)'(/,'..)%/./..&..$*,+-.-/*/++-%.+--+.),-,#,&(//.-,).-/./-.'./..///....-(*-/)/.#-$.//.%$*(*..())+/--/.../..(,+,..*).-/#.*'.+/-(/...&'+)/--.+'-))&$,'//*'/.-/(/.).,/,-,*.,/.*.".)..)./*/*/+/.'/%/&-),/..,)).(/--/&,*/..,*)-#+/.,-,./-**./%.//+/*)().(+-,,-.%'//../--.--/)*#..,*///.-   NM:i:58 ms:i:325    AS:i:325    nn:i:0  ts:A:+  tp:A:P  cm:i:65 s1:i:397    s

This error:
[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'

Means that you have not generated an index for a BAM file. This needs to be done for most downstream processing of BAM files:
# generate index
samtools index gill.sorted.bam
# count genes
htseq-count -f bam --strand=no gill.sorted.bam yyy.gtf > gill-counts.txt

Update (after second error):
I've got no idea what's going wrong; this is not an error that I would expect to see, and doesn't look related to PacBio input. I wonder if the error is related to the number of threads used in the sort command, given that it happens at record 79, and you've specified 80 threads. The next things that I would do would be to look at the input files around where the error is to try to work out what's going wrong:
samtools view gill.sorted.bam | head -n 85 | tail -n 10

Yes, i have used 80 threads. Do i need to reduce the thread number?

That depends on what machine you're running on. Attempting to run 80 threads on a machine that only has 4, for example, will slow down computation and may result in unexpected behaviour.


text:
Slow bwa shm "preloading" into memory
Coming from devops background, working on automating usual genomics pipelines, I am benchmarking existing pipelines in order to choose cloud resources properly.
Past week most of work was on testing bwa mem.
I've noticed that bwa takes about 10 minutes to start processing data. When monitoring bwa resource usage I've noticed that RAM is slowly building up to about 6GB. That is amount of reference indexed data.
After bwa loads the reference first time, any subsequent run is fast. But why this preloading is happening so slow? My reference data set is already indexed with bwa index.
When I monitor disk activity with iotop I noticed that disk read is only 5MB/s, when I run bwa mem alignment (only first time), which is very slow. Performance of the drive is lot more than that and a simple copy test confirms that.
Is it a "limitation" of bwa, is this slow preloading expected to be like that and I should just accept it and move on forward or there is something I don't understand and maybe someone can help me understand what could I do?
The same "slow preload" happens when I run bwa shm /path/to/ref.
Then bwa mem starts immediately. But again slow preload is something that is bothering me and I am trying to understand what can I do about it.

Issue was with EBS snapshot which holds reference indexed data. Issue is that when you run EC2 from snapshot, EBS storage performance is expected to be low until you optimize the volume. My assumptions and questions for two weeks, all "goes to the water".
Unfortunately I was not aware of this hence I asked the question out from ignorance.

Can you please add this as an actual answer and accept it? This will help the StackExchange algorithms to identify this as a question with and answer, and stop promoting it in the future to fish for more answers.

I will. Didn not want to do that imediatelly since it is a little bit awkward to answer own question = )

Yeah, I find it awkward as well, but you did already answer the question (as a comment to your own question). Might as well get a bit more credit for that. Also, it's worth thinking about comments as transient things that are there only to help people improve their work; they're not meant to be used for anything substantial that should hang around for a long time (beyond gratitude).

I'll add an answer to my own question, which was asked out of lack of experience with cloud (AWS) EC2 service.
The way I was testing bwa is from a custom AMI image created for our pipelines.
The image has "preloaded" reference data that we are using for alignement so that we spare some time downloading reference data from S3 and spare some S3 api calls since we have a lot of files used as reference data.
Simply our image is single 200GB EBS volume used with AWS ECS AMI where we stored our reference data and created a custom image for Batch out of it. Then we use that AMI for AWS Batch jobs.
Once I used iotop (after few days of testing) to monitor disk activity I've noticed that bwa loads index into memory very slow, cca 5MB/s. Syncing files from S3 was faster then using this snapshot approach.
For benchmarking I've used the same AMI manually launched. In the AWS documentation there is part which explains why and how you should optimize EBS when AMI is launched from a snapshot.
https://docs.aws.amazon.com/AWSEC2/latest/UserGuide/ebs-initialize.html
After spining up fresh EC2, syncing reference data and running bwa everything was clear. My EBS volume was slow and had to be initialized first.
That is all. My lack of EBS storage knowledge was the issue.
This helped me understand that I should not use custom AMI with ref data stored and will probably use EFS or simply syncing the data from S3 inside Docker entrypoint.

Answer from @user438383, converted from comment:
BWA has slow indexing due to the indexing size; this can be improved by using bwa-mem2:
https://github.com/bwa-mem2/bwa-mem2

We are happy to announce that the index size on disk is down by 8 times and in memory by 4 times due to moving to only one type of FM-index (2bit.64 instead of 2bit.64 and 8bit.32) and 8x compression of suffix array. For example, for human genome, index size on disk is down to ~10GB from ~80GB and memory footprint is down to ~10GB from ~40GB. There is a substantial reduction in index IO time due to the reduction and hardly any performance impact on read mapping.

Yes but it is not yet production ready. Talking from my very poor amount of reading on the bwa-mem2 topics.


text:
Using grep to get lines from a file (.tsv) that contain specific samples listed in a sample file (.tsv)
So I have a file, let's call it pcawg.tsv. It's formatted like the one down below.
Tumor sample ID  Normal sample ID  Mechanism
a                1                 deletion
b                2                 ecDNA
c                3                 ecDNA

I have a file of tumor sample ID's with the ecDNA mechanism called ecDNA_samples.tsv. It's formatted like this:
Tumor sample ID
a
b

To filter and extract the lines of pcawg.tsv that contain the tumor sample ID's in ecDNA_samples.tsv into a new file, I used this line of code below.
grep -f ecDNA_samples.tsv pcawg.tsv > ecDNA.vcf.tsv

Would this be the right code to use?

This has been asked on https://unix.stackexchange.com/questions/110645/select-lines-from-text-file-which-have-ids-listed-in-another-file

Is there a particular reason why you want to use grep to do this? It would be possible to make a more robust solution using R.

I kind of have to use this with Linux because I'm working with a Linux-based secure computing portal, since we have patient data which is confidential.

This may help:
cat ecDNA_samples.tsv | xargs -I % grep '^%\s' pcawg.tsv

The command above executes grep '^<line>\s' pcawg.tsv, for each line of ecDNA_samples.tsv. Which matches the lines taken from ecDNA_samples.tsv to the sequences of characters between line start and a white-space in pcawg.tsv.

So I ran this code as cat ecDNA_samples.tsv | xargs -I % grep '^%\s' pcawg.tsv>ecDNA_normal.tsv, but when I opened the file ecDNA_normal.tsv in my directory, the file was empty.

Use >> instead of > so an empty result at the end doesn't end up as the result instead of everything you've captured.

I see! Thank you!

@RamRS It worked.

Just in case an R solution is acceptable as well, this should be fairly robust code because it properly parses the TSV files and matches up column names:
#!/usr/bin/env Rscript

library(tidyverse)
read_tsv("pcawg.tsv") %>%
   inner_join(read_tsv("ecDNA_samples.tsv")) %>%
   write_tsv("filtered_pcawg.tsv")


text:
Seurat - single cell analysis
I have been doing single cell analysis based on seurat and I realized after each clustering and reclustering when I run SCTransform, in my code there is vars.to.regress = percent.mt and I was running this everytime. I removed MT- contaminated cells before.
My clustering seems fine though.
Do you think this is very critical thing and I have to analyze everything from scratch?

I think it would be a good idea to at least look at the distribution of mitochondrial gene fraction in your cells. If they mostly have low fractions (e.g. less than 15%), then I wouldn't expect the clustering to be influenced too much by them.
With SCTransform regressing on the percent.mt, I think it's okay to do post-transform filtering, but it'd be a good idea to exclude any cells that have high mitochondrial fraction prior to doing clustering (unless, for whatever reason, those are the cells you're interested in).


text:
calculating phylogenetic tree from pairwise distances
I have a file with a 3-column list of pairwise distances:
A B distance
A C distance
B C distance
...

And I would like to calculate a phylogenetic tree for all the nodes in the list. I've seen phylip neighbor and similar options, but it seems like the length of the ids is a limit, which I rather not have to deal with. What would be a straightforward way of getting a newick tree output out of an input as above?

I believe there are a number of ways to construct a tree metric from a distance metric. A very straightforward method, neighbor joining, is available in the Sciki-Bio package. A less straightforward option with more freedom is by using the scipy.cluster.hierarchy module to obtain a linkage matrix, then using the to_tree() method to obtain a scipy Tree object, and finally writing a script to convert a scipy Tree object to a newick string.
I've copied the example from scikit-bio below:
>>> from skbio import DistanceMatrix
>>> from skbio.tree import nj
>>>
>>> data = [[0,  5,  9,  9,  8],
...         [5,  0, 10, 10,  9],
...         [9, 10,  0,  8,  7],
...         [9, 10,  8,  0,  3],
...         [8,  9,  7,  3,  0]]
>>> ids = list('abcde')
>>>
>>> dm = DistanceMatrix(data, ids)
>>> # Or from a tsv file with header, no column 0
>>> reader = csv.reader(open(args.inputfile), delimiter="\t")
>>> x = list(reader)
>>> dm = DistanceMatrix(x[1:],x[0])
>>>
>>> newick_str = nj(dm, result_constructor=str)
>>> print(newick_str[:55], "...")
(d:2.000000, (c:4.000000, (b:3.000000, a:2.000000):3.00 ...

Note: Both methods require you to supply the input as a numpy distance matrix


text:
functions in seurat to calculate the gene count per cluster
I'm running FindAllMarkers
markers <- FindAllMarkers(scData_scTransform_PTX_only_res04_20PC_6000g_30mt_sb_removal_20210713_ann_exp0_counts, min.pct = 0.25, logfc.threshold = 0.25)

But was wondering if there is a function in seurat to get the raw (non-normalized) expression data from each cell type cluster

Hi, could you please add more context to this question so that we can understand why you want to find out this information? I've calculated cell counts per cluster, and visualised gene counts per cluster using scatter plots, but haven't yet run into a case where I'd need to work out gene count per cluster as a single statistic (whatever that means).

@mmpp could it be that you meant to compare expression profiles of some genes (by means of a boxplot, for instance), to evaluate the differences between two conditions, in the same cell type, i.e. a cluster?

are you trying to calculate pseudobulk average?


text:
Understanding the dot plot from seurat
I am working with single cell data and using seurat to analyze the results. Often in manuscript, we see the dotplots showing the expression of the marker genes or genes of interest across the different clusters. My question here is:
a. When comparing expression of a gene (lets say geneA) across cluster, for example, lets say for cluster 0, the average expression for geneA that is plotted, is it calculated considering all the cells in the sample or all the cells in the cluster0? If it is calculating the average expresison based on total number cells present in cluster 0 (here), is it correct to compare the expression of geneA across all the clusters?
Because among the different clusters, the number of cells are different, and the number of cells that are expressing geneA can vary alot across the clusters.
Thank you in advance!

Hi, could you please add more context to your question? An example image of a Seurat dotplot that demonstrates your issues would help a lot here.

Since it is showing mean expression among cells within a cluster, why would different number of cells across clusters worry you?


text:
General Question: Statistical Alternatives to Oncoprint
I am a computational genomics student who has recently joined a research lab as an RA. We often work with oncoprints (comuts/heatmaps), and we visualize molecular signatures by reorganizing the comut. Although this is the first-line approach, what are some statistical ways to identify molecular signatures? I have heard of NNF, clustering, etc, but I also heard a colleague mention using cox analysis (which I thought was for survival analysis?).
My Background:
While introductory statistics has been hammered into my head, I have less experience with higher level computational techniques, like PCA, clustering. I am familiar with the general purpose of the statistical approaches I've mentioned, but not confident in when or how to use them. For example, I know the general purpose of PCA, but I did not know that if you have categorical binary data and want to do unsupervised clustering, you should do logisticPCA. I would appreciate any advice, general or relevant to the question at hand. I am eager to be competent, confident, and contribute to my lab! Thank you :)

Can you tell us more about what you are trying to accomplish? Or link to examples? For example I see that comut and oncoprint are specific and distinct visualization softwares (https://github.com/vanallenlab/comut, https://bioc.ism.ac.jp/packages/3.2/bioc/vignettes/ComplexHeatmap/inst/doc/s8.oncoprint.html), Cox was a prominent statistician who worked on many techniques, not only survival analysis, but you seem to be using the words in very generic ways so I worry you mean something different. Can you be more specific?

I agree with Max, I have hard time understanding your question. You don't need to explain your presonal context, but rather clearly explain what exactly you do (perhaps with links), what you identify as a problem (why you can not continue doing so) and how do you imagine to solve it.


text:
Problem with Seurat reference mapping
I have 10X scRNA-seq multiome's 3' poly-A capture (scRNA-seq+ATAC-seq) from PBMC
Using Seuratreference mapping, I mapped my scRNA-seq part on reference PBMC (for cell type annotation)
I have also tried Seurat zimuth shiny app
But likely just a few of my cells has been mapped on the reference
Please look at this

Do you think this because my data is multiome or something else is going wrong?
Please any body has some experience help me
Thank you

Yes! Multiome seems to have cells with much higher mitochondrial gene content. I don't understand the biology behind this. But in using a cut off of 5% (i.e. remove cells with >5% MT genes) I lose over 80% of cells in multiome data. Also, here is the Signac tutorial where they use a 20% cut off for multiome data - https://satijalab.org/seurat/articles/weighted_nearest_neighbor_analysis.html ("WNN analysis of 10x Multiome, RNA + ATAC")

As I mentioned in my answers to similar questions, there are a lot of cells with very high mitochondrial gene content. It would be a good idea to work out why that is before doing any more downstream processing.
These cells will, of course, be filtered out by the thresholds you have set.


text:
How to read and edit a FASTA file with python using regular expression?
I have a file containing protein sequence information of 100 odd proteins. Each sequence starts with a header that looks like this:
>sp|Q9UDW1|QCR9_HUMAN Cytochrome b-c1 complex subunit 9 OS=Homo sapiens OX=9606 GN=UQCR10 PE=1 SV=3

I want to edit these header to look like this:
>Q9UDW1

I have tried the following code:
import re
with open ('WT.fasta', 'r+') as f:
    contents = f.readlines()
    for i in range(len(contents)-1):
        target = re.compile(r'>..\|(.*)\|.*', contents[i])
        target.sub(r'\1', contents[i])

I think I am not being able to call the right methods to do this. I have been trying to read the documentation: https://docs.python.org/3/library/re.html. I am not been able to understand it and feeling overwhelmed by the details. I would really appreciate if anyone could help or point me in the right direction. Thanks!

Do you absolutely need to use python? A simple sed would do this much faster.

Hi @RamRS, I am actually very new to bioinformatics. I still have not learnt about sed or awk commands :-)

You definitely should - these will save you a lot of time for tasks that are too small in get Python/R involved.

@RamRS, I would love learn about them. Would you have any suggestions about any course, or tutorials on these topics?

Just google around. The way to learn these would be using them on copies of files for tasks like your current question - start there and you'll find more and more places you can use them.

Thanks a lot :-)

In general, I usually recommend Biopython for any sort of fasta parsing, but the code below should work as well.
with open ('WT.fasta', 'r+') as in_f, open("WT_out.fasta", 'w') as out_f:
    for line in in_f:
        if line[0] == ">":
            out_f.write(">" + line.split('|')[1] + "\n")
        else:
            out_f.write(line)

Hi, The code executes, does not show any error messages. But the file looks the same as before!

@Noob which file are you looking at? The code creates a new file WT_out.fasta

Oh sorry! My bad :/... I am an absolute idiot!

You could do this with awk by processing the header:
$ awk -v FS="|" '{ if ($0 ~ /^>/) { printf(">%s\n", $2) } else print $0 }' in.fa > out.fa

For example:
% echo -e ">sp|Q9UDW1|QCR9_HUMAN Cytochrome b-c1 complex subunit 9 OS=Homo sapiens OX=9606 GN=UQCR10 PE=1 SV=3\nACCGT" | awk -v FS="|" '{ if ($0 ~ /^>/) { printf(">%s\n", $2) } else print $0 }'
>Q9UDW1
ACCGT

Hi @AlexReynolds, Thanks a lot! I will try it.


text:
DESEQ2 Question about results()
I am currently learning to perform Differential Analysis via DESEQ2 R Package, and I believe I've made progress, able to format the data correctly [maybe] for DDS(). When I run the results function to see the output, the data seems fine, as can be seen below:
log2 fold change (MLE): Group RA vs C  Wald test p-value: Group RA vs
C  DataFrame with 6 rows and 6 columns
       baseMean log2FoldChange     lfcSE      stat    pvalue      padj
      <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric> Gnai3   95.1978      0.3458366  0.275366  1.255916  0.209146  0.999383
Cdc45   21.9942     -0.0545177  0.346199 -0.157475  0.874870  0.999383
H19     29.0133     -1.4954220  1.516390 -0.986173  0.324048  0.999383
Narf    16.6967      0.2571059  0.370880  0.693232  0.488164  0.999383
Cav2    63.6356      0.5060019  0.354969  1.425480  0.154018  0.999383
Klf6    65.8469      0.2851175  0.265755  1.072858  0.283335  0.999383

However, when I run summary() or any other function, they do not seem to behave correctly, not counting any stats. Same with MAplot not coloring genes that should be marked red:
out of 13184 with nonzero total read count adjusted p-value < 0.1 LFC
> 0 (up)       : 0, 0% LFC < 0 (down)     : 0, 0% outliers [1]       : 150, 1.1% low counts [2]     : 0, 0% (mean count < 1)

I'm confused as to why this is occurring and was wondering if anyone had ideas.
EDIT: Added More Code for context
sergey <- (read.csv("cleanmatrix25+.csv", sep = "\t", header = T,
    check.names = F, row.names = "Gene_names"))

meta <-(read.csv("metadata25+.csv", sep = "\t", row.names = 1))

dds <- DESeqDataSetFromMatrix(countData = sergey, colData = meta,
    design = ~Group)

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

dds$Group <- relevel(dds$Group, ref = "C")

dds <- DESeq(dds)

res <- results(dds)

head(res)

summary(res)

plotMA(res, ylim=c(-2,2))

EDIT2:
I sorted result by p-value, so all of these should be red on the following MAplot due to having a p-value < 0.1

You do not have any significantly differential genes after FDR correction according to padj, simple as that, which means that either there are none or your sample size is too small (underpowered) given the observed variability of your data. As the MA-plot shows some logFCs that could qualify as large enough to be called DEG it probably comes down to a too small sample size and quite some variability between the replicates of the groups. A PCA can help diagnosing whether there is simply lots of per-group variability or other factors such as batch effects that need to be taken care of. Please see the DESeq2 vignette for code examples or Bioc packages such as PCAtools for this.
Edit: plotMA by the way colors based on padj, not p-value alone.

Thank you so much! I will try the PCA to see if there are any significances


text:
How to check if a fasta file and a GTF file fit and form a valid pair?
I'd like to know if there is a simple way to check the concordance between a GTF file and a fasta file (that I use as mapping reference) ?
It may be a dumb question, but I suspect discrepancies between my 2 files and I don't want to check this manually

Answer to this specific question --> turn back to the person who provided me the GTF and complain ! :P

Answer from @ATpoint, converted from comment:
Check whether chromosome identifiers are the same and both files contain the same chromosomes. You can also check that the end coordinate of the GTF never exceeds the length of the respective chromosome. Maybe you can also check whether the first three nucleotides of every CDS indeed start with a start codon. If this is all true then I guess I would trust this pair.

You might use agat_sq_filter_feature_from_fasta.pl from AGAT
 DESCRIPTION
 This script is a kind of annotation filter by sequence name. It goes through the gff annotation features and remove those
 that are not linked to a sequence from the fasta file provided. The
 match between sequence name in the fasta file and the 1st column of
 the gff3 file is case sensitive.


text:
Adjusted P-value from DESeq2
What does the padj really mean? I believe it's a BH method, where it attempts to quantify the likelihood of a "false positive" but I am not sure. Also, what does it mean if the resulting calculations nearly all have the same padj values despite having different p-values?

It would help to see the data themselves (p-values and adjusted p-values), or at least a general idea of how many hypotheses and the part of the p-value spectrum. For high p-values FDRs ("adjusted p-values", https://en.wikipedia.org/wiki/False_discovery_rate) will be very similar as they will mostly indicate support for the null hypothesis, i.e. very high false discovery rates. I would suggest reading about DESeq2 stats (section "Multiple test correction"): https://hbctraining.github.io/DGE_workshop/lessons/05_DGE_DESeq2_analysis2.html

If we look at documentation of DESEQ2 and search "adjusted p-value", we find the section "Multiple test correction". In this section they discuss Benjamini-Hochberg ("BH") false discovery rate (FDR) correction procedures:

FDR/Benjamini-Hochberg: Benjamini and Hochberg (1995) defined the concept of FDR and created an algorithm to control the expected FDR below a specified level given a list of independent p-values. An interpretation of the BH method for controlling the FDR is implemented in DESeq2 in which we rank the genes by p-value, then multiply each ranked p-value by m/rank.

Earlier, they write:

m = total number of tests

They further add:

So what does FDR < 0.05 mean? By setting the FDR cutoff to < 0.05, we’re saying that the proportion of false positives we expect amongst our differentially expressed genes is 5%. For example, if you call 500 genes as differentially expressed with an FDR cutoff of 0.05, you expect 25 of them to be false positives.

I would suggest reading Benjamini and Hochberg (1995) to learn more details if desired.


text:
PyMol: Output Polar Contacts Between Chains to Text File
I am new to PyMOL but have a very specific task that I need to do.
I have a PDB structure file of protein homo-oligomer, and I want to use PyMol to determine polar contacts between chain A and all other chains.
However, I not only want to visualize these polar contacts (which is easy to do via PyMOL GUI), but I also want to be able to store/save pairs of interacting residues in a dictionary or whatever may be the most convenient data structure to do this. The end goal to create a sort of "interaction map" in a csv or text file where each interacting pair of residues has its own line.
I would ideally like to do this via a Python script where I call the PyMol API or through a PyMol script - whichever is easier.
I have seen a similar question on Stack Overflow but found it unhelpful because I not only want to select and visualize these polar contacts but also want to save/record them in a list or dictionary (and then to a text or csv file), which the above question does not address.

Please do not shotgun post to multiple sites. (Post on biostars: https://www.biostars.org/p/9484676/)

You should give a try to ChimeraX's hbonds or contacts.

So you can install PyMOL as a standalone, but you can also install it as a bona fide Python 3 module via conda:
conda install -c schrodinger -c conda-forge pymol-bundle

It can be installed in other ways —apt-get or brew or even compiled, but the latter is excruciatingly painful.
In your python notebook you can do:
import pymol2
with pymol2.PyMOL() as pymol:
    pymol.cmd.load('structure_file.pdb', 'foo') # file
    # pymol.cmd.fetch('PDB_CODE') # PDB code
    # pymol.cmd.read_pdbstr(foo_str, 'foo') # string
    ...

where the ellipsis can be any API cmd.
This context manager approach is to avoid the problem of namespace pollution as the session (e.g. loading structures without removing the old ones etc.) is the same but you don't see what is happening. Also PyMOL() is parallel the parallel mode (so mutagenesis wizard does not work), while in singleton mode (pymol2.SingletonPyMOL), you get issues if you try to re-initialise etc.
Anyway. To "select" all interacting atoms at the interface you can do:
n = pymol.cmd.select('chain A around 3')

Expand includes the selection. Around does not.
This adds a selection to the session. So what we actually want is a pythonic instance to play with. So:
atoms = pymol.cmd.get_model('chain A around 3')
where atoms is a list of chempy.Atom objects. With attributes such as resi, resn, chain etc. Do dir on an atom for more...
One attribute is coord, which contains a length 3 list of floats.
NB. The above selection gives not a full residue, whereas 'byres chain A around 3' does. So if you want to make a list of all interacting residues:
interactors = [{'chain': atom.chain, 'index': atom.resi, 'number': atom.resn} for atom in pymol.cmd.get_model('name CA and byres chain A around 3')]

An iteration can give the atom to atom distances.
Do note that this is all proximity, so you will have to screen out hydrogens and carbons to get actual polar contacts.

After some searching, I discovered the command 'cmd.find_pairs(selection1, selection2, mode, angle, distance)'. If you pass it two selections, it will return pairs of atoms within the given criteria. This seems to be better than using "around" because it takes into account atom orientation - which is important when inferring hydrogen bonds.


text:
Limma decideTests function: what kind of multiple hypothesis testing correction does parameter "method" involve?
What kind of multiple hypothesis testing correction does method="global" do in Limma's decideTests function? According to the documentation:

method="global" will treat the entire matrix of t-statistics as a
single vector of unrelated tests.

and

method="global" is useful when it is important that the same
t-statistic cutoff should correspond to statistical significance for
all the contrasts.

What does this really mean? Does it use a Bonferroni correction?

More context from the docs, looking at the arguments to decideTests():

method: character string specifying how genes and contrasts are to be combined in the
multiple testing scheme. Choices are "separate", "global", "hierarchical"
or "nestedF".

adjust.method: character string specifying p-value adjustment method. Possible values are "none",
"BH", "fdr" (equivalent to "BH"), "BY" and "holm". See p.adjust for details.

These are two different options.
One option (adjust.method) chooses the multiple test correction, for which the default is BH/FDR.
The other option (method) decides how the multiple test correction is applied, for example can it be applied separately to different genes or across contrasts. method = global in this context just means that every single test statistic is considered at the same time. An alternative to this would be to consider a nested design.
More intuitively, method = separate applies the correction to each column of your matrix independently, such that in a $m x n$ matrix of statistics, the "number of tests" is $m*n$ for global and $m$ for separate.
For understanding the effect of the number of tests on FDR, I suggest consulting the FDR wiki page. Basically, using a method other than global will increase power.

What do m and n correspond to here?

m would be number of rows of the test statistic matrix, n number of columns. So the total number of test stats in a single column is m, the total number of test stats in the full matrix is $m * n$. Columns of matrix correspond to slices of the analysis as defined by limma, IIRC. So you can choose to correct every possible test at once, or you can correct by slice of the analysis.


text:
Subsetting a Seurat object based on colnames
I'm trying to subset my seurat object based on colnames. I have gone ahead and labeled each cluster and now I want to subset all the colnames that are in Cancer_human for human_colnames and all the barcodes that are not in Cancer_human for mouse_colnames but get an error.
human_colnames = colnames(scData[,scData$cell_ann == "Cancer_human"])
mouse_colnames = colnames(scData[,!scData$cell_ann == "Cancer_human"])

Error in `[.Seurat`(scData, , scData$cell_ann == "Cancer_human") : Incorrect number of logical values provided to subset cells

I'm filtering on colnames not rownames. I want all the barcodes from the Cancer_human.

The code could only make sense if the data is a square, equal number of rows and columns. At the moment you are getting index from row comparison, then using that index to subset columns. Also, please provide a reproducible example data for testing, dput(myData).

@zx8754 i agree for the reproducible example (consider the pbmc data set that comes with Seurat) but not the square argument (it would be correct for matrices) but Seurat objects contain multiple things and is implemented such that for a Seurat object so, so[, i] and so[[j]][i] both give information about cell i.

Did you tried the subset function?
# assuming that scData@meta.data looks like this
#         orig.ident  ...  cell_ann
# ATTCAA  sample1          Cancer_human
# ATTTAT  sample1          Cancer_human
# ...
# AAATCA  sampleX          Cancer_mouse
# ATTTAA  sampleX          Cancer_mouse

Idents(scData) <- "cell_ann"
human_subset <- subset(scData, idents = "Cancer_human")
mouse_subset <- subset(scData, idents = "Cancer_mouse")
mouse_subset <- subset(scData, idents = "Cancer_human", invert = TRUE) # alternative way

P.S. It is not clear what you are trying to do. You should add to your post how your data looks like.

If you want the column names to be preserved, you can subset the column names as a vector:
human_colnames = colnames(scData)[scData$cell_ann == "Cancer_human"]
mouse_colnames = colnames(scData)[!scData$cell_ann == "Cancer_human"]


text:
Splitting .bam files into separate samples for tophat2
I have RNA seq data of 12 samples (paired end, so 24 files)  and I ran tophat2 for alignment on all of them which gave me one bam file for all the samples. Now I need separate bam files for all the samples. Can I spilt them from the one I have or I have to do tophat2 for all the samples separately? If yes, how to do it.
I need them for FRASER (splicing detction) I would also appreciate if anyone used FRASER before could help me a bit to get started. I am working with plant data for the first time.
The command I used
tophat2 -p 40 \
    -o tophat_output_final bowtie_indexes/Medicago1/Medicago \
    NRNSA_S103_R1_001.fastq.gz,NRNSB_S107_R1_001.fastq.gz,\
      NRNSC_S111_R1_001.fastq.gz,NRSA_S104_R1_001.fastq.gz,\
      NRSB_S108_R1_001.fastq.gz,NRSC_S112_R1_001.fastq.gz,\
      RNSA_S105_R1_001.fastq.gz,RNSB_S109_R1_001.fastq.gz,\
      RNSC_S113_R1_001.fastq.gz,SRAB_S106_R1_001.fastq.gz,\
      SRB_S110_R1_001.fastq.gz,SRC_S114_R1_001.fastq.gz \
    NRNSA_S103_R2_001.fastq.gz,NRNSB_S107_R2_001.fastq.gz,\
      NRNSC_S111_R2_001.fastq.gz,NRSA_S104_R2_001.fastq.gz,\
      NRSB_S108_R2_001.fastq.gz,NRSC_S112_R2_001.fastq.gz,\
      RNSA_S105_R2_001.fastq.gz,RNSB_S109_R2_001.fastq.gz,\
      RNSC_S113_R2_001.fastq.gz,SRAB_S106_R2_001.fastq.gz,\
      SRB_S110_R2_001.fastq.gz,SRC_S114_R2_001.fastq.gz

Did you assign read groups while aligning? Can you add the command that you used with tophat2?

Hi Wouter De Coster, I have added the command in the question itself. I don't know what you mean by assigning read groups?

The way you aligned them you can not split them in separate samples. You should have run Tophat2 for each sample separately.

Thank you, Wouter De Coster I have 3 replicates of each treatment. so 4 treatments x 3 replicates =12 files.
Should I merge the replicates or just go for separate alignment for each sample?

Yes, separate alignments for each (biological) replicate


text:
Transcriptomics vs Proteomics for understanding molecular responses
I am planning some experiments to investigate how mono cultures of microbes respond to different nutrient concentrations. I want to use an -omics based approach to understand how, on a molecular level, the microbes utilise these nutrients by looking for up regulation of associated genes. With the idea looking out for upregulation in mind, I was going to go for transcriptomics. But now I am beginning to realise that regulatory changes in genes does not necessarily translate to the molecular utilisation of a nutrient by a cell. To that end, perhaps I would be better using quantitative proteomics to capture changes in the concentration of functional proteins which are actually interacting with the nutrients. However, I worry that by relying on proteomics I might miss important regulatory changes for microbe species where the proteome is less well defined.
I should say that it's an either/or situation, there is only budget for one.
I'd appreciate any thoughts on the best approach.
Thanks.

You might reach a greater audience at biostars.org. Do you have access or a potential collaboration with a group that can run the high-throughput proteomics? RNA-seq is simple, one might even say trivial these days, both if you do the libraries yourself or just send the RNA to a facility while proteomics is a bit less standard I think, at least we had get in contact with a collaborator to get some proteome-wide experiments going, both because we lack experience and most importantly the machinery (HPLC, MS and everything around it).

My gut feeling is that you will be figuring out a lot more causality using the transcriptomic approach, but that might be just because I am a lot more familiar with that type of data. This is a really open-ended question I think, a lot more suited to discussion forums than Q&A sites.

I agree that this is a discussion question. Proteomics is less available but also less sensitive, albeit requiring less corrections analysis-wise. For transcription factors and other low copy number protein transcriptomics is the best option. Bacteria do not go around sequestering mRNA as much as Eukaryotes, which is nice, but some phyla, produce compounds that make them problematic for RNA work (e.g. the osmolyte in Firmicutes spores (something-picolinate)), so proteomics > transcriptomics for Geobacillus etc. Also... there is metabolomics if enzymes are involved (but way less sensitive).

Also, before starting any -omics, check the quality of your reference genome if it is a non-standard species. This sounds unimportant, but if your have an auto-annotated reference genome, you will spend more time fixing it than anything else —e.g. pathway X is missing enzyme Y, but in the operon with WYZ there is a mystery enzyme, which guilty by association would suggest is a candidate, but you have to get more evidence.

@MatteoFerla Thanks for this. When you say quality, are you referring to coverage? The species I'm interested in has a reference genome with 39% coverage, is that likely to be sufficient?

No, I meant the quality of the annotations: the data is meaningless if too many of the genes are genes of unknown function or worse not what the automatic annotation claimed. I am assuming you have full coverage. reference genome with 39% coverage is unusable.


text:
Why is cutadapt discarding some paired end reads?
So, I am using cutadapt to remove two primers from some paired end reads.
I have 2 files for each sample, i.e. sample_R1.fastq and sample_R2.fastq.
I have primer_1 and primer_2 and both of them can be found in sample_R1.fastq and sample_R2.fastq. In other words, if I search for primer_1 and primer_2 in sample_R1.fastq, the sum of the matches is the same as the number of the reads in sample_R1.fastq. Same for R2.
To properly trim them what I do is the following, through R:
setOne <- paste("-g", primer_1, "-a", primer_1)
setTwo <- paste("-G", primer_2, "-A", primer_2)
system2(cutadapt, args = c(setOne, setTwo, "-n", 2,
                           "-o", outOne, "-p", outTwo,
                            inputOne, inputTwo, "--minimum-length", 1))

The command works fine except for one thing. On the one hand cutadapt should remove those sequences which, after the trimming, will have a length < 1, as specified by the command. However, some reads are removed even though they are longer than 1 and I would like to know why this happens, since I can't find a specific reference mentioning that cutadapt removes sequences.
If I am getting this straight, cutadapt only removes the primer sequence from the beginning of the read, then removes the equivalent characters in the quality field of the read that was trimmed and, finally, it returns the trimmed read with the quality field adjusted to the new length. If the length is lower than the threshold, the read is discarded. If it is discarded in R1, it will be discarded in R2 and viceversa.
I give an example of read that is removed, this one from R1:
@some:sample:number 1:N:0:1
ACCGTCACGTCATCAGCAGCATCGACTCCTGTTTGATCCCCGTGCCTGTACTCAAGTCAACCAGTACACGCTTTTCAAGCTGCCTTCTGCCGTATTTCACACTCAGCGCATTTGCAATTGGTGTTCTTTGCACCACATATCTATCAGTTAAGCCACCGCTACTCCGCCTACCTCGTACTGACTTAATCGACCGCCTACGCACCCTTTAAACCCAATTCGAGACCGTCAATTCAATCATACTTTAGTAAGTG
+
BBBBBBBBBBBBGGGGHEGHFHA@EGGGHHHHHHHHHGHHHHHHHHGGGGGHGHGHGGGHHHHHGDGGGFHGGGGHGGGGGGHGHGFHHHHHHHHHHGGGHHHGHGHHHHHHHGFFGGGGGGGGGGGGGGGGGFFBDFDFFFFFFFBFFFFFFFFFHGHGGHHHHHEGGGGGHHHHGHHHHHGHHHHHHHHHHHHGGGHGHHGGGGEFGGHHHHHHHHHHHHHGHHHHHHHHHHHFHFGHHHHGGGGGHHG

this one from R2:
@some:sample:number 2:N:0:1
CTGCTGCTATCACTAGCTAAAGCTTTAACGTAGGGCAATTGACGAAGGTGAAATGCATAGATATGTGATGACGGTGTAGCGTTGGCGCTGGGGTTTAAAGGGTGCGTAGTGACTTCACAGAACACCAATTGGCAGCTTACTAAAGCATGAGTACAGACGAGGTAGGCGGAATTGCGGTCGATTAAGTCAGTTGTGAAATACGGCAGCTTAACTGTTGGCAAGCGTTGTCCGGCACTGGTTTTATTCGAGGT
+
AAAAA@BFFFFFGFGCGGFDDC/CDHF2FGFHHG//<<CCC1FGHH1<1>11C<FF/0BFFFFFFF/:DEEEFFHGG0AEG/B9B./-.@9GHGEFC>EE<FHEHEGHDHGGFHGGBGGGHGGDGH0CDDGFFHHCC--.:CGCFG??GGAEGGHFHHBFFFEFB./B..;/;:/999=AAA../BB/;//////9/EFFCBDAA..990ADD=99/;/BFGGGGGGGGGFHHHHHHGHDGHHFFDHHHHG

and these are the primers that I am removing:

ACCGTCACGTCATCAGCAGCATCGAC
CTGCTGCTATCACTAGCTAAAGCTTTA

matching the reads I provided. Clearly, if I trim them manually, the reads length is way greater than 1 so, why are these reads removed from the final cutadapt output file?
EDIT:
I might have solved the issue, but I asked on the cutadapt github issues page, just to make sure.


text:
How to convert multiple single-end bam files to fastq using samtools
Hi I am trying to convert bam files generated from Ion Torrent Proton sequencing to fastq format so that I can upload them to KBase for analysis.
The files are named as follows:
01_thu159_IX_RNA_001.bam 01_thu163_IX_RNA_001.bam 01_thu164_IX_RNA_001.bam 01_thu190_IX_RNA_001.bam
So I tried:
samtools fastq *.bam > *.fastq
However, that did not work. But, I know it is possible to use a placeholder for the name in code by using something like $f or I have also seen %N%. I am just not sure how to basically use those placeholders to be able to quickly apply a for loop, which will convert the files to fastq format without having to do each one individually. I have 28 samples each containing 4 bam files. Any advice on how I can do this would be greatly appreciated.

I am not familiar with samtools however your example, which I believe is run through the terminal, is printing all the output into a single file and not into many, as you may think. I guess that by using > *.fastq you meant to create many fastq files, which is not what you're doing.
also I found this link that may help you https://onestopdataanalysis.com/bam-to-fastq/ and here it looks like they process one file at a time. if this is the case you just need a for loop to run the command for each file

Hi there, yes I was running that command in the terminal. So the idea basically was I don't want to have to type the name of the individual files in each case. For example I would like to know how to write the for loop so that it converts the bam file to fastq format and maintains the same file name just the format is changed. Otherwise I will have to type out each name individually as such: samtools fastq 01_thu159_IX_RNA_001.bam > 01_thu159_IX_RNA_001.fastq and this is tedious considering how many files I have.

ok, so your question actually was about how to run a command on multiple files at once, am I right? If yes, then here you can find some answers to that specific question:
https://unix.stackexchange.com/questions/80851/how-do-i-run-a-command-on-multiple-files
https://unix.stackexchange.com/questions/80851/how-do-i-run-a-command-on-multiple-files

Yes, that is correct. Thank you, I truly appreciate the assistance.

of course, if you are not able to write the necessary code, you can still ask a question regarding that specific point. However, I believe there should be enough information in the two answers I linked to help you.

Using GNU parallel:
ls *.bam | parallel -j 2 "samtools fastq {} > {.}.fastq"

Use -j to control the number of parallel jobs. That only works for single-end data though as paired-end data must first be sorted or collated by name.

Thank you very much, I really appreciate it. My data is single-end reads, as it was generated using Ion Torrent, so this should work well.

@Justin1609 Yes, I guess that should work with the above command without problems.


text:
Where to find publicly available Sanger chromatography data?
I am looking for any public available databases where I can download Sanger chromatography data ideally in Ab1 or SCF file format. I need extremely large amounts of this, since I want to use it for machine learning.

@ATpoint , can you link me here any ML based solution for basecalling if you know any ? I found only 2-3 papers in the literature

Why do you want to reinvent the wheel on this? This is an old technology that doesn't require any new tools, not least of which because its use is decreasing every year.

@DevonRyan I can fully unterstand what do mean but Sanger DNA sequencing is not going away any time soon, also we are looking for inhouse Basecalling method/software which will be independent of any licencing

I must agree with the others. Furthermore, I'd say writing a script that can deal with frameshifted sequences in the PCR pool (the major issue with basecalling) is easier, but that is application specific as in general Sanger is used for validation and the data has to be perfect. E.g. I wrote a tool to call secondary peaks to test the diversity of error prone PCR pools, which is a super specific application.

There are some interesting trends in a Sanger trace, i.e. peak height distribution, which tell you about how well the experiment was (concentrations, gunk etc.), but that again is not too interesting to be heavily used —unless say implemented in Benchling or other major tool, i.e. the tricky part of a new tool is getting to be adopted.
Furthermore, about licensing the Ab1 files are what is spat out by the machine and that is where the licensing/compression comes into play, not in basecalling as one can very easily do the basecalling from the raw data within them ignoring the calls.

NCBI Trace Archive
The NCBI Trace Archive is a permanent repository of DNA sequence chromatograms (traces), base calls, and quality estimates for single-pass reads from various large-scale sequencing projects. Contaning more than 2 billion traces.


text:
Snakemake rule issue with bcftools view filtering
I am creating a bioinformatics pipeline using Snakemake doing everying from adapter trimming of raw sequencing data in .fq to mapping, SNP calling and SNP filtering.
The issue i am having is when i try to filter specific chromosomes of the bcf files, using the --regions. The rule i have created below works and does create a filtered bcf file.
rule bcftools_filter:
    input:
        bcf = "/dir/{sample}_aligned_call.bcf",
    output:
        "/dir/{sample}_aligned_call_filter.bcf"
    threads: 20
    shell:
        "bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} -o {output}"

However i would like to add the --regions parameter to the bcftools command. I have tried the shell part of the rule to several things
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions <chr1,chr2> -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {chr1,chr2} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {{chr1,chr2}} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {{chr1} {chr2}}
 {{chr2}} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions ['chr1','chr2'] -o {output}

But it doesn't work, and i am having trouble finding any suggestions to this specific issue.
I know it works if i just run bcftools filter in my command line, so using --regions should work.
Could anyone have any suggestions how to add --regions to the snakemake rule?

Have you tried expand() function? https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html?highlight=expand#the-expand-function

No i havent tried that, i didn't realize you could potentially use the expand() function for the shell part. But i will look into it

Something like the following would work:
rule bcftools_filter:
    input:
        bcf = "/dir/{sample}_aligned_call.bcf",
    output:
        "/dir/{sample}_aligned_call_filter.bcf"
    params:
        regions="chr1,chr2"
    threads: 20
    shell:
        "bcftools view -Ob --threads {threads} -m2 -M2 --regions {params.regions} -v snps {input.bcf} -o {output}"

You could also include the regions as part of the output file name, rather than specifying it as a param.

Would a simple --regions \"chr1,chr2\" in the shell command not work?

You don't need to even quote it. However having it as a param allows you to set it elsewhere (e.g., the command line).


text:
Alignment statistics

Given a multi-fasta multiple sequence alignment file, is there a quick way/tool to calculate average pairwise sequence identity?
Even better, given a whole genome alignment file (.maf or even .hal) is there a quick way/tool to calculate average pairwise sequence identity across the individual local multiple sequence alignments?

Thanks for reminding me, and apologies for not accepting earlier. I found what I was looking for in maftools, but agree the answer below was also useful.

There are tons of methods that do something similar to this. Is it necessary that it be alignment rather than a heuristic such as k-mer hashing?
If k-mer hashing is ok, then I would suggest mash or sourmash. There are additionally newer fancy techniques.
For just an old-school program that computes pairwise identities, PHYLIP has programs DNADIST and PROTDIST that exactly compute identities. Similar questions on biostars suggest that PHYLIP's utilities are recommended even relatively recently elsewhere. But it looks like Clustal Omega also computes identities for free.
For microbial genomes a common statistic is average nucleotide identity (ANI), which may be what you are looking for. There are a variety of tools to do this, here is one, and here is a webserver that will do it for you.


text:
Assumptions of batch effect removal
What does removing a batch effect (e.g. with limma::removeBatchEffect) assume about the batch effect?
Does it assume simply a constant batch effect for each level?
E.g. if the levels are
Date1
Date2
Date3

Then one constant will be subtracted from all samples on Date3, another constant from all samples of Date2 and so on?
Or is it more complicated?

You're basically subtracting a constant per-gene per-level. The relevant portion of the code is:
fit <- lmFit(x,cbind(design,X.batch),...)
beta <- fit$coefficients[,-(1:ncol(design)),drop=FALSE]
beta[is.na(beta)] <- 0
as.matrix(x) - beta %*% t(X.batch)

where x is your input matrix, design is the design matrix and X.batch is the matrix of batch and covariates. So you're fitting your data with the combined batch/design matrix (line 1), extracting out the batch coefficients (line 2), and subtracting out the expected effect (line 4). This is simple enough and remarkably effective for its intended use. You can get a bit fancier by using PCA or something else to instead allow individual samples to be more/less affected by a batch-effect, but added complexity has its own issues.


text:
Where to find rat genome for STAR alignment?
I'm trying to find the rat genome to carry out STAR alignment for several FASTQ files.  Where would I find these files?

For example from Ensembl: https://www.ensembl.org/Rattus_norvegicus/Info/Index
Can find other links via a search engine of choice.


text:
Build protein-protein interaction prediction deep learning model
I'm an undergraduate biology student and my thesis is on designing a deep learning architecture to predict whether two proteins interact or not given their primary sequences. I have read some papers with different approaches to the problem including the paper [https://academic.oup.com/bioinformatics/article/35/14/i305/5529260] uses PIPR model which is the state of the art model for this problem in 2019 and recently in 2021 the paper [https://www.biorxiv.org/content/10.1101/2021.01.22.427866v1] uses D-SCRIPT which surpass PIPR with the aspect of the small dataset and has better generalization feature for many model organisms.
For my problem, I have two questions

can you suggest to me the approach to build a model for my thesis. Should I follow these papers to build a model, how my model will be different from the models in the papers?
Which architecture or model recently that the literature considers the most proper for the protein-protein interaction prediction approach?

I very much appreciate any help from everyone!

Albeit to make a full structure... Check out https://github.com/sokrypton/ColabFold it's a mod of the AlphaFold2 notebook. Albeit not mentioned by AlphaFold2, the Rosettafold paper talks about complexes. It's not benchmarked. If two proteins don't interact the replicates will differ.

I’m voting to close this question because it is opinion-based and not on-topic. Thesis guidance should be done by a local supervisor, not by strangers on the internet. This has also been asked over at biostars.org (but I cannot find the link) where an experienced member added a good response poiting out that the aim of this endeavour seems to go way beyond of what an undergraduate can/should do.

Here is the post mentioned by @ATPoint: https://www.biostars.org/p/9485815/


text:
Extract read names and the associated nucleotides on specific positions from a BAM file (in R)
Let's assume I have a BAM file and several positions that I would like to examine more closely in this alignment. My goal is to find out whether these positions are on the same reads and which nucleotides each read has in these two positions.
So my output should be a table of all reads that occur in the positions and list the associated nucleotides.
The output should look like this (with 980, 1350 and 1400 the first three positions; I made up the read names):

read name
980
1350
1400
...

ReAdN4m1L03
A
T

...

ReAdN4m2L02
A
T
G
...

ReAdN4m3L01

T
G
...

...
...
...
...
...

At the end I would like to have a table that links the nucleotides via the reads.
So far, I have tried to read in the Bam file at the desired positions:
#Positions of interest
pos <- c(980, 1350, 1400)

#path to BAM file
path <- pathToBamFile

#name of reference used for mapping
refname <- "myreference"

#loading the bam file
bamFile <- BamFile(path)
gr <- GRanges(seqnames = refname, ranges = IRanges(start = pos, end = pos))
params <- ScanBamParam(which = gr, what = scanBamWhat())
aln <- scanBam(bamFile, param = params)

With the following code I can access the reads, but how do I get the corresponding nucleotides in the individual positions?
pos.no <- c(1:length(pos))

i <- 1
aln[[pos.no[i]]]$qname

I have already tried to work with the CIGAR strings aln[[i]]$cigar with the sequences aln[[i]]$seq in order to extract to the nucleotides manually, but that led to more and more problems.
Is there any easy way to extract the nucleotides?
I have searched for a long time now on the internet for an easy solution but without success. Hopefully someone will help me here! Many thanks for your time in advance.

Do you have to use R? This is more convenient if you use the pileup interface from pysam (perhaps Rsamtools provides a similar functionality).

@DevonRyan but does pileup give you an output where the nucleotides are linked to the reads?

Yes, it does, at least in non-R interfaces (i.e., python and C).

This answer does it in Python, perhaps there are some ideas there: https://bioinformatics.stackexchange.com/questions/7401/access-base-aligned-to-particular-reference-position


text:
Error when using awk command to trim sequence files
I am attempting to edit my fastq sequence files with an awk command from a published article. The sequence files are all contained in one directory where I am running the awk command from the terminal in a Linux environment. There are 4 fastq files in the directory
01_thu159_IX_RNA_001.fastq
01_thu163_IX_RNA_001.fastq
01_thu164_IX_RNA_001.fastq
01_thu190_IX_RNA_001.fastq

I would like to trim all these files at once with the awk command so I wrote this for loop:
for filename in ./*.fastq*
do
awk -v s=10 -v e= 40 ' {if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;} ' ./${filename} >  ./${filename}_trimmed
done

However, I am getting the following error:
awk: fatal: cannot open file ` {if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;} ' for reading (No such file or directory)

I then thought maybe it was an error with the for loop so I tried the command as per the instructions in the article
awk -v s=10 -v e= 40 '{if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;}' ./01_thu159_IX_RNA_001.fastq > ./Trim_01_thu159_IX_RNA_001.fastq

But this also produced the same error:
awk: fatal: cannot open file `{if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;}' for reading (No such file or directory)

I have checked the syntax for the awk command and it should be correct. Could anyone please give me advice as to what I am doing wrong and how to fix this issue?

You have a whitespace here:
e= 40

awk thinks that 40 is your rule and the rest of the line is a filename.

Thank you very much really appreciate it, I didn't notice that. Will give it a try now again.


text:
What does the following diagram represent?
Suppose, the following is a diagram of a protein's polypeptide chain:

What does this diagram represent?

What are the letters A, E, M, W, L, N, S, etc. represent?  (I suppose these are amino acid codes. Right?)
Which one is the c-alpha atom here?
What are the sequence numbers i-2, i-1, i, i+1, i+2, etc. represent? Why there are negative and positive increments from i?

The following is a related diagram:

Yes, 2. it's a schematic with a whole residue was a circle —from Voet and Voet IIRC. 3. Residue i is the one in focus, the reference. i-1 is the preceding residue. Have a gander at https://en.m.wikipedia.org/wiki/Turn_(biochemistry) where this type of relative notation is important.

is this from a school or university assignment?

@user3479780, this is not an assignment. this is a personal study for data analysis.

Yes, the letters are amino acids
It's a schematic with a whole residue as a circle —possibly from Voet and Voet.
Residue i is the one in focus, the reference. i-1 is the preceding residue. Have a gander at en.wikipedia.org/wiki/Turn_(biochemistry) where this type of relative notation is important.

The diagram shows how secondary structure (helices, sheets etc.) affect the distance between the residues. In a helix a residue 3 positions down will be close, but not in a sheet.
However, for assignment of secondary structure it may be easier to look at the φ and ψ angles:

α helix: phi=-57.8, psi=-47.0
π helix: phi=-57.1, psi=-69.7
3.10 helix: phi=-74.0, psi=-4.0
β sheet: phi=-139, psi=+135

Although these have a range and you can get π bulges within α helices —important for conformation switching for pore opening (capsaicin binding to TRPV1 is a classic example).


text:
How likely is it to find primers sequences in already-trimmed reads?
So, I am analysing 18S amplicon data (fastq files) to be able to eventually investigate the taxa composition.
After removing primers (using cutadapt) from both R1 and R2 (assuming primers are found at the beginning of my sequences), I am still able to detect, apparently at random, a few reads that contain my primers (i.e. a part of the sequence matching the primer) in, again apparently random, positions.
I am able to detect almost all my primers (16 out of 24, both forward and reverse) in at least one sequence and I am wondering how common (or how weird) this is. I guess it is important to consider that my primers contains a few, i.e. from 1 to 3, wild-card bases that may match multiple other bases (apparently called whobble bases), which makes it easier to match a random sequence (there are multiple possibilities for an actual match).
Problem is that the presence of primers in the sequences may affect downstream analyses (e.g. with DADA2 they are very clear about the fact that primers have to be removed) in several ways and I was wondering if these sequences I am finding may be interfering somehow of if I can live with it.

The primers will be at the 3'-end of the reads. Make sure you check for them there.

@DevonRyan yes, that I know and that I did so no primers at 3'. I guess it just happens to find them inside the reads, as random matches, but I wasn't expecting that many actual matches. I guess it's statistics, then.

Interesting, have you tried to map a read with a sequence matching a primer to the reference, or even simply blasting it? Does the reference contain the "primer sequence"? ("" as we don't know if it's a homologous sequence or not)

@KamilSJaron apparently this is something that should not be happening and I'm trying to understand why it does.


text:
Encoding different DNAs in pySBOL3?
I want to encode two different DNAs using pySBOL3. First dna is an engineered region usually composed by a transcriptional unit (TU) or a set of TUs. The second is vector which is a basically a plasmid that contains a set of the previously described dna. This is for differentiate between engineered on different plasmids.
Now for dna I used as type: sbol3.SBO_DNA and role: sbol3.SO_ENGINEERED_REGION and for vector I used as type: sbol3.SBO_DNA and role: SO:0000155 (plasmid).
I would like to have feedback on this choice. For example the plasmid could also be an engineered plasmid SO:0000637 or circular plasmid SO:0002211 or just circular SO:0000988.
Here is my actual implementation:
for dna_id in params.get('dna'):
    #this needs to be handlded by vector
    dnas = Dna.object().filter(dna_id)
    dna = sbol3.Component(dnas.name(), sbol3.SBO_DNA)
    dna.roles = sbol3.SO_ENGINEERED_REGION
    dna.name = dnas.name
    dna.description = dnas.description
    if dnas == None:
        print('No dnas specified')

for vector_id in params.get('vector'):
    vectors =Vectors.object().filter(vector_id)
    vector = sbol3.Component(vectors.name(), sbol3.SBO_DNA)
    vector.roles = SO:0000155
    vector.name =vectors.name
    vector.description = vectors.description()

Your choices match the same choices that I've made when designated types for inserts and vectors.
The more specific plasmid types are reasonable to use as well, but choosing just "plasmid" is certainly safe. That's why I have preferred to use that in my code because (1) my code might sometimes be applied to a different sort of plasmid, (2) I don't currently have a purpose for making the finer distinction, and (3) circular DNA would be indicated by SO.circular (SO:0000988) in any case.
Side note: I would suggest using the tyto package so that your code can use names like SO.plasmid in your code instead of opaque symbols like SO:0000155.


text:
Differential miRNA expression using RPM
I have microRNA (miRNA) expression data in RPM. I would like to do differential gene expression between two groups.
How can I do this?
I have considered edgeR and DESeq2 in R, but it looks like they require counts and I cannot use RPM in these. Or can I convert RPM to counts?
Can anyone help?

This has been asked many times before. Here is an answer from the edgeR/limma author. Approach 2 in the toplevel question is probably the best choice here if you don't have counts: https://support.bioconductor.org/p/125144/#125156 so an analysis with limma starting from the RPMs on log scale.

Thank you. I will try that.


text:
How to determine NCBI's SRS google cloud bucket or AWS bucket
It appears that AWS and GCP host SRA data and it's beneficial to grab the data from that source when running on GCP for example.
Given an SRR accession like SRR1929796 https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1929796
Is there a metadata API / source that can tell me which google cloud bucket to retrieve the data from like they show on the website?
gs://sra-pub-crun-4/SRR1929796/SRR1929796.1

It also sounds like fastq-dump and fasterq-dump can be configured with Google Cloud credentials to download that for me.
I'm just wondering is there a way to determine the google cloud storage bucket path it's in gs://sra-pub-crun-4?

It sounds like it to me too. I have never directly accessed NCBIs stuff via Google buckets and I don't know what permissions they are setting for direct access.
So ... I briefly checked it out
gsutil ls gs://sra-pub-crun-4

and 'yes' its possible but ...

BadRequestException: 400 Bucket is a requester pays bucket but no user project provided.

'Requester pay' luckily I wasn't logged in :-)
What the message implies given access to your (my) GCP account and they'll bill you (me). I don't think it will be much, particularly for a 'list' request, however I ain't gonna find out for a data set I don't need. I do doubt I'm going to get billed for 'list' request, however I've never paid 'list' request I prefer to keep it that way. In contrast, NCBI's front-end is free and thats what they want, because they could of course set the bucket so they pay.
So in answer the question, if you wanna pay go ahead ... dial up a GCP account link it to pay.

Here's my first answer
Please note I'm not really certain what information the OP seeks here, e.g. gs:// locations type stuff. If the question is where's the SRS database on every gs:// path ... I don't know.
I think the OP wants an API ... Google Cloud certainly use APIs by the bucket load (no pun intended but I thought it was funny). I don't use them because they are excessively verbose and I prefer SDK and gcloud tools, of which gsutil is one.

It also sounds like fastq-dump and fasterq-dump can be configured with Google Cloud credentials to download that for me.

If I have understood correctly and that is all you want then sure
gsutil cp gs://sra-pub-crun-4/SRR1929796/SRR1929796.1 .

will do it, if permitted (or in hindsight prepared to pay) and there will be an API version of this.

Background
Installation via Google's Cloud SDK, via  https://cloud.google.com/storage/docs/gsutil_install#install