text:
Converting VCF format to text for use with PLINK and understanding column mapping
I successfully completed Nature PRS tutorial, which is based on PLINK. Turning to my real data, I downloaded ukb-d-20544_1.vcf.gz. Now I'm facing the problem that I seem to be unable to use it in PLINK or find the correct data format to download at all, and I am a bit confused about the correct general summary statistics-download procedure.
1. Convert .vcf format to a .txt file, PLINK can understand
In the tutorial they use a .txt.gz base data file aka summary statistics that looks like this:
CHR BP  SNP A1  A2  N   SE  P   OR  INFO    MAF
1   756604  rs3131962   A   G   388028  0.00301666  0.483171    0.997886915712657   0.890557941364774   0.369389592764921
1   768448  rs12562034  A   G   388028  0.00329472  0.834808    1.00068731609353    0.895893511351165   0.336845754096289
...

Having this .vcf.gz (please see excerpt at the very bottom), I unzipped and, in R, used MungeSumstats::format_sumstats("ukb-d-20544_1.vcf", ref_genome = "GRCh37")[1], to get something out of it similar to what they use in the .txt file in the tutorial. However, the column names won't match exactly. I know that I can use BETA instead of OR, and FRQ might correspond to MAF, but I can't seem to find the imputation information score INFO. Here is what MungeSumstats::format_sumstats produced:
SNP CHR BP  A1  A2  END FILTER  FRQ BETA    SE  LP  N   P
rs12238997  1   693731  A   G   693731  PASS    0.116734    0.00058531  0.000756397 0.357493    117706  0.439042942120051
rs371890604 1   707522  G   C   707522  PASS    0.0983981   0.000793177 0.000848559 0.456023    117706  0.349926634613021
...

I need to have the columns corresponding to MAF and INFO as per the tutorial.
The next step in the tutorial involves filtering based on specific criteria, such as MAF and INFO values. However, I can't proceed because these columns are missing.
I would like to primarily use the terminal and only turn to R for statistical analysis. Could you guide me on the correct approach to convert the .vcf file to .txt so that I can continue following the tutorial with PLINK?
2. General approach
In addition to addressing the conversion issue, I'm also seeking guidance on the appropriate data sources for genetic analysis. I'm curious about whether OpenGWAS is a reliable source for publication or if there are other recommended sources like dbGaP. Additionally, what is the ideal data format to obtain?

The beginning of the received .vcf file looks like this:
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##FORMAT=<ID=ES,Number=A,Type=Float,Description="Effect size estimate relative to the alternative allele">
##FORMAT=<ID=SE,Number=A,Type=Float,Description="Standard error of effect size estimate">
##FORMAT=<ID=LP,Number=A,Type=Float,Description="-log10 p-value for effect estimate">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Alternate allele frequency in the association study">
##FORMAT=<ID=SS,Number=A,Type=Float,Description="Sample size used to estimate genetic effect">
##FORMAT=<ID=EZ,Number=A,Type=Float,Description="Z-score provided if it was used to derive the EFFECT and SE fields">
##FORMAT=<ID=SI,Number=A,Type=Float,Description="Accuracy score of summary data imputation">
##FORMAT=<ID=NC,Number=A,Type=Float,Description="Number of cases used to estimate genetic effect">
##FORMAT=<ID=ID,Number=1,Type=String,Description="Study variant identifier">
##META=<ID=TotalVariants,Number=1,Type=Integer,Description="Total number of variants in input">
##META=<ID=VariantsNotRead,Number=1,Type=Integer,Description="Number of variants that could not be read">
##META=<ID=HarmonisedVariants,Number=1,Type=Integer,Description="Total number of harmonised variants">
##META=<ID=VariantsNotHarmonised,Number=1,Type=Integer,Description="Total number of variants that could not be harmonised">
##META=<ID=SwitchedAlleles,Number=1,Type=Integer,Description="Total number of variants strand switched">
##META=<ID=TotalControls,Number=1,Type=Integer,Description="Total number of controls in the association study">
##META=<ID=TotalCases,Number=1,Type=Integer,Description="Total number of cases in the association study">
##META=<ID=StudyType,Number=1,Type=String,Description="Type of GWAS study [Continuous or CaseControl]">
##SAMPLE=<ID=ukb-d-20544_1,TotalVariants=9938650,VariantsNotRead=0,HarmonisedVariants=9938650,VariantsNotHarmonised=0,SwitchedAlleles=209,StudyType=Continuous>
##contig=<ID=1,length=249250621,assembly=HG19/GRCh37>
##contig=<ID=2,length=243199373,assembly=HG19/GRCh37>
...
##contig=<ID=22,length=51304566,assembly=HG19/GRCh37>
##contig=<ID=X,length=155270560,assembly=HG19/GRCh37>
##contig=<ID=Y,length=59373566,assembly=HG19/GRCh37>
##contig=<ID=MT,length=16569,assembly=HG19/GRCh37>
##contig=<ID=GL000207.1,length=4262,assembly=HG19/GRCh37>
##contig=<ID=GL000226.1,length=15008,assembly=HG19/GRCh37>
...
##contig=<ID=GL000192.1,length=547496,assembly=HG19/GRCh37>
##gwas_harmonisation_command=--json /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1_data.json --ref /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/QC/genomes/b37/human_g1k_v37.fasta; 1.1.1
##file_date=2019-11-25T15:12:05.258952
##bcftools_annotateVersion=1.9-74-g6af271c+htslib-1.9-64-g226b4a8
##bcftools_annotateCommand=annotate -a /mnt/storage/home/gh13047/mr-eve/vcf-reference-datasets/dbsnp/dbsnp.v153.b37.vcf.gz -c ID -o /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1.vcf.gz -O z /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1_data.vcf.gz; Date=Mon Nov 25 15:37:00 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  ukb-d-20544_1
1       692794  rs530212009     CA      C       .       PASS    AF=0.111693     ES:SE:LP:AF:SS:ID       0.000379398:0.000798869:0.197332:0.111693:117706:1_692794_CA_C
1       693731  rs12238997      A       G       .       PASS    AF=0.116734     ES:SE:LP:AF:SS:ID       0.00058531:0.000756397:0.357493:0.116734:117706:rs12238997

For the first part of question 2:  if there are other recommended sources like dbGaP.
Please refer to this Databases of GWAS summary statistics:

Database
Content

GWAS Catalog
GWAS summary statistics and GWAS lead SNPs reported in GWAS papers

GeneAtlas
UK Biobank GWAS summary statistics

Pan UKBB
UK Biobank GWAS summary statistics

GWAS Atlas
Collection of publicly available GWAS summary statistics with follow-up in silico analysis

FinnGen results
GWAS summary statistics released from FinnGen, a project that collected biological samples from many sources in Finland

Pheweb.jp
GWAS summary statistics of Biobank Japan and cross-population meta-analyses

https://www.nature.com/articles/s43586-021-00056-9/tables/3

Thanks for pointing to that paper, very helpful +1

However, the column names won't match
exactly. I know that I can use BETA instead of OR, and FRQ might
correspond to MAF, but I can't seem to find the imputation information
score INFO.

Yes. There is not a standard format that everyone reports GWAS summary statistics in, unfortunately. First, check the paper that you are sourcing it from to see if they already filtered on imputation quality. Do ctrl+f and look in the methods section for imputation information (you may also try searching for imputation tools, like minimac or Beagle).
The good news is that if the sumstats have many alleles, imputation accuracy doesn't matter much. Errors get averaged out and each error only has a small effect since there are many alleles. Usually it is expected that reported sumstats will be pre-QC'd in some way, especially if they lack INFO. It would be safe to proceed without the INFO filter step.

The next step in the tutorial involves filtering based on specific
criteria, such as MAF and INFO values. However, I can't proceed
because these columns are missing.

Sumstat formatting is whatever the author wants, and can be a bit annoying to parse. However, FRQ is basically MAF.
See here (ctrl+f "frq"):
https://bioconductor.riken.jp/packages/3.14/bioc/vignettes/MungeSumstats/inst/doc/MungeSumstats.html
https://academic.oup.com/bioinformatics/article/37/23/4593/6380562

I would like to primarily use the terminal and only turn to R for
statistical analysis. Could you guide me on the correct approach to
convert the .vcf file to .txt so that I can continue following the
tutorial with PLINK?

Did you try unzipping and simply replacing .vcf with .txt in the filename?

General approach In addition to addressing the conversion issue, I'm also seeking guidance on the appropriate data sources for genetic
analysis. I'm curious about whether OpenGWAS is a reliable source for
publication or if there are other recommended sources like dbGaP.
Additionally, what is the ideal data format to obtain?

This aligns with your first question. If you use a database, sumstats will be standardized (usually) to the database's common format, so there will be less confusion over inconsistent column names. The GWAS catalog pretty good in my opinion, and all of them are acceptable.


text:
missed cleavages combinations generation
I am looking for an R based solution to the following problem:
Trypsin does not cleave perfectly in proteomics experiments. So starting from a protein sequence FASTA I want to generate all possible peptides which COULD be generated if trypsin was allowed n number of "missed cleavages"
In informatics terms I would describe this as skipping a split during a strsplit() function
Example to clarify:
We start out with the protein sequence: GRGKA which is usually split by strsplit() (digested by trypsin) after every K and R
"GRGKA" %>% strsplit(., "(?<=K|R)", perl = TRUE)

[[1]]
[1] "GR" "GK" "A"

now imagine that trypsin did not cleave perfectly, for example missing always one K or R
"GRGK" "GKA"

would be generated in this case. Or if trypsin digested this protein really badly and would end up missing two cleavages we would end up with the sequence:
"GRGKA"

I would love to have a generalizable solution that allows to generate all combinations of sequences that arise with n number of missed cleavages.
If you need any more information please feel free to ask, sorry if I missed something!
jay.sf's requested example:
"FLGERTCTRNERQP"

perfect digest:
"FLGER" "TCTR" "NER" "QP"

one missed cleavage all combinations:
"FLGERTCTR" "TCTRNER" "NERQP"

two missed cleavages all combinations:
"FLGERTCTRNER" "QP" "FLGER" "TCTRNERQP"

three missed cleavages (there are a total of three K and Rs so trypsin would miss all)
"FLGERTCTRNERQP"

For longer sequences I think this becomes a frameshift like problem: If we count the K and Rs in the sequence:
FLGERTCTRNERQPFLGERTCTRNERQP
FLGE1TCT2NE3QPFLGE4TCT5NE6QP
FLGER_TCTR_NER_QPFLGER_TCTR_NER_QP (perfect digest)
FLGERTCTR_NERQPFLGER_TCTRNER_QP (one missed cleavage and always the uneven R and Ks)
FLGER_TCTRNER_QPFLGERTCTR_NERQP  (one missed cleavage and always the even R and Ks)
So I think for 1 missed cleavage we have two frames over the whole sequence and for two missed cleavages we would have three frames and so on....

Say "FLGERTCTRNERQP" would be a candidate, it has three possible splits. Can you show all possible results for 1, 2, and 3 skipped splits?

FLGERTCTRNERQP
perfect digest:
FLGER TCTR NER QP
one missed cleavage all combinations:
FLGERTCTR TCTRNER NERQP
two missed cleavages all combinations:
FLGERTCTRNER QP FLGER TCTRNERQP
three missed cleavages (there are a total of three K and Rs so trypsin would miss all)
FLGERTCTRNERQP

@KlemensFröhlich the answer below looks cool, you've enough rep to vote - definitely recommend it.

Essentially we want to combine 1, 2, ..., n - 1 splits, and if the variance of split lengths is greater than one we have a second solution. Here is a function that,  after the strsplit, first creates indices cb to later paste the respective splits together. I'm not quite sure what the maximum number of splits is, but this might be a start.
mccg <- function(x) {
  ## missed cleavages combinations generation ++++++++++++++++++++++++++++++++++
  gs <- \(n, k, length=n*k, right=FALSE) {
    ## helper, similar to `gl()` but w/o factorization
    if (right) {
      cbr <- gs(n, k, n, right=FALSE)
      rev(max(cbr) - cbr + 1L)
    } else {
      rep_len(rep.int(seq_len(n), rep.int(k, n)), length)
    }
  }
  ## split string after [K, R]
  spl_x <- el(strsplit(x, '(?<=K|R)', perl=TRUE))
  len <- length(spl_x)
  res <- if (len <= 1L) {
    ## if [K, R] at the end or missing
    x
  } else {
    ## nested list of combinations
    cb <- lapply(1:(len - 1L) |> {\(.) setNames(., . - 1L)}(), \(i, .len=len) {
      cb <- list(gs(.len, i, .len))
      if (i != 1L && var(tabulate(el(cb)))) {
        ## second solution, starting from right
        c(cb, list(gs(.len, i, .len, right=TRUE)))
      } else {
        cb
      }
    }) |> c(setNames(list(rep_len(1L, len)), len - 1L))
    ## paste according to sequences
    rapply(cb, \(x) vapply(unname(split(spl_x, x)), paste, collapse='',
                           character(1)),
           how='list')
  }
  return(res)
}

Usage
Two splits:
mccg(x='GRGKA')
# $`0`
# [1] "GR" "GK" "A"
#
# $`1`
# $`1`[[1]]
# [1] "GRGK" "A"
#
# $`1`[[2]]
# [1] "GR"  "GKA"
#
#
# $`2`
# [1] "GRGKA"

Three splits:
mccg(x='FLGERTCTRNERQP')
# $`0`
# [1] "FLGER" "TCTR"  "NER"   "QP"
#
# $`1`
# [1] "FLGERTCTR" "NERQP"
#
# $`2`
# $`2`[[1]]
# [1] "FLGERTCTRNER" "QP"
#
# $`2`[[2]]
# [1] "FLGER"     "TCTRNERQP"
#
#
# $`3`
# [1] "FLGERTCTRNERQP"

More splits:
mccg(x='FLGERTCTRNERQPFLGERTCTRNERQP')
# $`0`
# [1] "FLGER"   "TCTR"    "NER"     "QPFLGER" "TCTR"    "NER"     "QP"
#
# $`1`
# $`1`[[1]]
# [1] "FLGERTCTR"  "NERQPFLGER" "TCTRNER"    "QP"
#
# $`1`[[2]]
# [1] "FLGER"       "TCTRNER"     "QPFLGERTCTR" "NERQP"
#
#
# $`2`
# $`2`[[1]]
# [1] "FLGERTCTRNER"   "QPFLGERTCTRNER" "QP"
#
# $`2`[[2]]
# [1] "FLGER"          "TCTRNERQPFLGER" "TCTRNERQP"
#
#
# $`3`
# $`3`[[1]]
# [1] "FLGERTCTRNERQPFLGER" "TCTRNERQP"
#
# $`3`[[2]]
# [1] "FLGERTCTRNER"     "QPFLGERTCTRNERQP"
#
#
# $`4`
# $`4`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTR" "NERQP"
#
# $`4`[[2]]
# [1] "FLGERTCTR"           "NERQPFLGERTCTRNERQP"
#
#
# $`5`
# $`5`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTRNER" "QP"
#
# $`5`[[2]]
# [1] "FLGER"                   "TCTRNERQPFLGERTCTRNERQP"
#
#
# $`6`
# [1] "FLGERTCTRNERQPFLGERTCTRNERQP"

Simulation
Here are some toy data,
set.seed(42)
n <- 1e3
min_chr <- 5; max_chr <- 50
strings <- replicate(n, {
  repeat {
    smp <- sample(setdiff(LETTERS, c("B", "Z", "J", "O", "U", "X")),
                  sample(min_chr:max_chr, 1L), replace=TRUE)
    if (any(grepl('K|R', smp))) break
  }
  paste(smp, collapse='')
}) |> setNames(nm=_)

over which the execution of the function can be tested in a loop.
r <- lapply(strings, \(x) tryCatch(mccg(x), error=\(e) NA))

hi jay
Thank you very much for looking into this! This already looks promising, but the main goal would not be to get a function that can go as deep as possible but rather that can be limited at any level that the user wants to choose.
So basically it would be ideal to have a function that allows to set the missed cleavages to 0, 1, 2, etc and always receive all the combinations up to this point. Do you think your function could be modified to control the depth of allowed missed cleavages? Best Klemens

source: Stack Exchange

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

sort: votes


text:
Generating gene signature for sample classification and survival analysis patient cohort
This is from this paper
The LSC17 score is calculated for each patient as a linear combination of GE of these 17 genes weighted by regression coefficients that were estimated from the training data as follows: LSC17 score = (DNMT3B × 0.0874) + (ZBTB46 × −0.0347) + (NYNRIN × 0.00865) + (ARHGAP22 × −0.0138) + (LAPTM4B × 0.00582) + (MMRN1 × 0.0258) + (DPYSL3 × 0.0284) + (KIAA0125 × 0.0196) + (CDK6 × −0.0704) + (CPXM1 × −0.0258) + (SOCS2 × 0.0271) + (SMIM24 × −0.0226) + (EMP1 × 0.0146) + (NGFRAP1 × 0.0465) + (CD34 × 0.0338) + (AKR1C3 × −0.0402) + (GPR56 × 0.0501). As above- and below-median scores in the training cohort were associated with adverse and favourable cytogenetic risk, respectively, a median threshold was used to discretize scores into high and low groups.
I'm unable to understand or apply the median threshold to create patient group specifically the R implementation here.
In order to address the problem my small dummy data-frame
This is my dummy data
set.seed(123)
nr1 = 20; nr2 = 8; nr3 = 6; nr = nr1 + nr2 + nr3
nc1 = 16; nc2 = 8; nc3 = 10; nc = nc1 + nc2 + nc3
mat = cbind(rbind(matrix(rnorm(nr1*nc1, mean = 1,   sd = 0.5), nr = nr1),
                  matrix(rnorm(nr2*nc1, mean = 0,   sd = 0.5), nr = nr2),
                  matrix(rnorm(nr3*nc1, mean = 0,   sd = 0.5), nr = nr3)),
            rbind(matrix(rnorm(nr1*nc2, mean = 0,   sd = 0.5), nr = nr1),
                  matrix(rnorm(nr2*nc2, mean = 1,   sd = 0.5), nr = nr2),
                  matrix(rnorm(nr3*nc2, mean = 0,   sd = 0.5), nr = nr3)),
            rbind(matrix(rnorm(nr1*nc3, mean = 0.5, sd = 0.5), nr = nr1),
                  matrix(rnorm(nr2*nc3, mean = 0.5, sd = 0.5), nr = nr2),
                  matrix(rnorm(nr3*nc3, mean = 1,   sd = 0.5), nr = nr3))
)
mat = mat[sample(nr, nr), sample(nc, nc)] # random shuffle rows and columns
rownames(mat) = paste0("gene", seq_len(nr))
colnames(mat) = paste0("LSC", seq_len(nc))

They have generated the LSC17 score by doing lasso,here in this case will create a dummy 5 gene score such as this
LSC5 score = (gene1 × 0.0874) + (gene6 × −0.0347) + (gene10 × 0.00865) + (gene22 × −0.0138) + (gene28× 0.00582)

Now to generate or categorize the samples into above or below the median score.
How to implement the above step in R which Im unable to do it?
One suggestion I received is this
 combined into a single variable when you make LSC17. After that, it's just quantile(LSC17,probs=c(0.3333,0.6666))

Which I unsuccessfully followed as I don't know how to implement it in the whole data frame.
Any help or suggestion would be really appreciated
I got an elegant solution which I asked here If there are other ways in R to do this I would love to know

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


text:
shiny app for single cell cluster export issue
I'm trying to create a shiny app for single cell data visualization and then label cells with genes till this step it works, the final step when I try to export cluster it is not exporting the desired cluster or i might be doing something wrong which Im not sure
This the code
library(Seurat)
library(BiocParallel)
library(reticulate)
library(Matrix)
library(shiny)
library(dplyr)

# Set the maximum file upload size to 3 GB
options(shiny.maxRequestSize = 3e+9)

# Define the UI
ui <- fluidPage(
  titlePanel("Single-cell RNA-seq Analysis"),
  sidebarLayout(
    sidebarPanel(
      fileInput("h5ad_file", "Upload H5AD File", accept = ".h5ad"),
      actionButton("run_button", "Run Analysis"),
      selectInput("gene_input", "Select Gene for Feature Plot:", choices = NULL),
      downloadButton("export_button", "Export Cluster")
    ),
    mainPanel(
      plotOutput("feature_plot")
    )
  )
)

# Define the server
server <- function(input, output, session) {

  # Function to perform the analysis
  performAnalysis <- function(h5ad_file) {
    sc <- import("scanpy")
    adata <- sc$read_h5ad(h5ad_file)

    counts <- t(adata$layers["counts"])
    colnames(counts) <- adata$obs_names$to_list()
    rownames(counts) <- adata$var_names$to_list()
    counts <- Matrix::Matrix(as.matrix(counts), sparse = TRUE)

    data <- t(adata$layers["winsorized"])
    colnames(data) <- adata$obs_names$to_list()
    rownames(data) <- adata$var_names$to_list()
    data <- Matrix::Matrix(as.matrix(data), sparse = TRUE)

    seurat <- CreateSeuratObject(counts)
    seurat <- SetAssayData(seurat, "data", data)
    seurat <- AddMetaData(seurat, adata$obs)

    seurat[["percent.mt"]] <- PercentageFeatureSet(seurat, pattern = "^MT[-\\.]")
    seurat <- NormalizeData(seurat)
    seurat <- FindVariableFeatures(seurat, nfeatures = 500)
    seurat <- ScaleData(seurat)
    seurat <- ScaleData(seurat, vars.to.regress = c("nFeature_RNA", "percent.mt"))
    seurat <- RunPCA(seurat, npcs = 50)
    seurat <- RunUMAP(seurat, dims = 1:20)
    seurat <- FindNeighbors(seurat, dims = 1:20)
    seurat <- FindClusters(seurat, resolution = 1)
    cl_markers <- FindAllMarkers(seurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = log(1.2))
    top_gene <- cl_markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)

    return(seurat)
  }

  # Reactive value to store the Seurat object
  seurat_object <- reactiveVal(NULL)

  # Perform analysis when the Run Analysis button is clicked
  observeEvent(input$run_button, {
    req(input$h5ad_file)
    h5ad_file <- input$h5ad_file$datapath
    seurat <- performAnalysis(h5ad_file)
    seurat_object(seurat)

    # Update gene choices for feature plot
    gene_choices <- rownames(seurat@assays$RNA@counts)
    updateSelectInput(session, "gene_input", choices = gene_choices)
  })

  # Render the feature plot
  output$feature_plot <- renderPlot({
    seurat <- seurat_object()
    if (!is.null(seurat)) {
      gene <- input$gene_input
      FeaturePlot(seurat, features = gene, label = TRUE)
    }
  })

  # Export cluster when the Export Cluster button is clicked
  observeEvent(input$export_button, {
    seurat <- seurat_object()
    if (!is.null(seurat)) {
      cluster_data <- cbind(Cluster = seurat@meta.data$cluster, seurat@assays$RNA@counts)
      write.csv(cluster_data, "cluster.csv", row.names = TRUE)
      updateActionButton(session, "export_button", label = "Download Cluster", href = "cluster.csv")
    }
  })

  # Function to end the session
  endSession <- function() {
    stopApp()  # This will end the Shiny session and shut down the server
  }

  # Call the endSession function when the session ends (browser window closed)
  session$onSessionEnded(endSession)
}

# Run the Shiny app
shinyApp(ui = ui, server = server)

The seurat object which is generated its structure is like this
str(seurat)
Formal class 'Seurat' [package "SeuratObject"] with 13 slots
  ..@ assays      :List of 1
  .. ..$ RNA:Formal class 'Assay' [package "SeuratObject"] with 8 slots
  .. .. .. ..@ counts       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. ..@ i       : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ...
  .. .. .. .. .. ..@ p       : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
  .. .. .. .. .. ..@ Dim     : int [1:2] 13714 2638
  .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
  .. .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..@ x       : num [1:2238732] 1 1 2 1 1 1 1 41 1 1 ...
  .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ data         :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. ..@ i       : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ...
  .. .. .. .. .. ..@ p       : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
  .. .. .. .. .. ..@ Dim     : int [1:2] 13714 2638
  .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
  .. .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..@ x       : num [1:2238732] 1.64 1.64 2.23 1.64 1.64 ...
  .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ scale.data   : num [1:500, 1:2638] 1.5265 -0.4982 -0.0324 -0.3129 -0.4631 ...
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:500] "MRPL20" "RER1" "TNFRSF9" "SRM" ...
  .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. ..@ key          : chr "rna_"
  .. .. .. ..@ assay.orig   : NULL
  .. .. .. ..@ var.features : chr [1:500] "PPBP" "LYZ" "S100A9" "IGLL5" ...
  .. .. .. ..@ meta.features:'data.frame':  13714 obs. of  5 variables:
  .. .. .. .. ..$ vst.mean                 : num [1:13714] 0.00341 0.00114 0.0019 0.00114 0.00682 ...
  .. .. .. .. ..$ vst.variance             : num [1:13714] 0.0034 0.00114 0.00189 0.00114 0.00678 ...
  .. .. .. .. ..$ vst.variance.expected    : num [1:13714] 0.00365 0.00114 0.00197 0.00114 0.00748 ...
  .. .. .. .. ..$ vst.variance.standardized: num [1:13714] 0.933 0.992 0.963 0.992 0.906 ...
  .. .. .. .. ..$ vst.variable             : logi [1:13714] FALSE FALSE FALSE FALSE FALSE FALSE ...
  .. .. .. ..@ misc         : list()
  ..@ meta.data   :'data.frame':    2638 obs. of  10 variables:
  .. ..$ orig.ident     : Factor w/ 1 level "SeuratProject": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ nCount_RNA     : num [1:2638] 2419 4903 3147 2639 980 ...
  .. ..$ nFeature_RNA   : int [1:2638] 779 1352 1129 960 521 781 782 790 532 550 ...
  .. ..$ n_genes        : num [1:2638] 781 1352 1131 960 522 ...
  .. ..$ percent_mito   : num [1:2638] 0.0302 0.0379 0.0089 0.0174 0.0122 ...
  .. ..$ n_counts       : num [1:2638] 2419 4903 3147 2639 980 ...
  .. ..$ louvain        : Factor w/ 8 levels "CD4 T cells",..: 1 3 1 2 5 4 4 4 1 6 ...
  .. ..$ percent.mt     : num [1:2638] 3.02 3.79 0.89 1.74 1.22 ...
  .. ..$ RNA_snn_res.1  : Factor w/ 8 levels "0","1","2","3",..: 1 3 1 6 4 1 5 5 5 6 ...
  .. ..$ seurat_clusters: Factor w/ 8 levels "0","1","2","3",..: 1 3 1 6 4 1 5 5 5 6 ...
  ..@ active.assay: chr "RNA"
  ..@ active.ident: Factor w/ 8 levels "0","1","2","3",..: 1 3 1 6 4 1 5 5 5 6 ...
  .. ..- attr(*, "names")= chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  ..@ graphs      :List of 2
  .. ..$ RNA_nn :Formal class 'Graph' [package "SeuratObject"] with 7 slots
  .. .. .. ..@ assay.used: chr "RNA"
  .. .. .. ..@ i         : int [1:52760] 0 102 187 274 457 833 1147 1503 1506 1550 ...
  .. .. .. ..@ p         : int [1:2639] 0 17 24 30 33 40 61 83 85 96 ...
  .. .. .. ..@ Dim       : int [1:2] 2638 2638
  .. .. .. ..@ Dimnames  :List of 2
  .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. ..@ x         : num [1:52760] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. ..@ factors   : list()
  .. ..$ RNA_snn:Formal class 'Graph' [package "SeuratObject"] with 7 slots
  .. .. .. ..@ assay.used: chr "RNA"
  .. .. .. ..@ i         : int [1:245992] 0 69 102 172 187 189 248 274 276 315 ...
  .. .. .. ..@ p         : int [1:2639] 0 79 118 179 253 323 424 541 574 650 ...
  .. .. .. ..@ Dim       : int [1:2] 2638 2638
  .. .. .. ..@ Dimnames  :List of 2
  .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. ..@ x         : num [1:245992] 1 0.0811 0.2903 0.1765 0.1765 ...
  .. .. .. ..@ factors   : list()
  ..@ neighbors   : list()
  ..@ reductions  :List of 2
  .. ..$ pca :Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
  .. .. .. ..@ cell.embeddings           : num [1:2638, 1:50] -3.544 -5.892 -3.544 10.139 0.932 ...
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..$ : chr [1:50] "PC_1" "PC_2" "PC_3" "PC_4" ...
  .. .. .. ..@ feature.loadings          : num [1:500, 1:50] 0.02701 0.16566 0.17077 -0.00786 -0.03674 ...
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:500] "PPBP" "LYZ" "S100A9" "IGLL5" ...
  .. .. .. .. .. ..$ : chr [1:50] "PC_1" "PC_2" "PC_3" "PC_4" ...
  .. .. .. ..@ feature.loadings.projected: num[0 , 0 ]
  .. .. .. ..@ assay.used                : chr "RNA"
  .. .. .. ..@ global                    : logi FALSE
  .. .. .. ..@ stdev                     : num [1:50] 5 3.39 3.2 2.78 1.94 ...
  .. .. .. ..@ key                       : chr "PC_"
  .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
  .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ]
  .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ]
  .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ]
  .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ]
  .. .. .. ..@ misc                      :List of 1
  .. .. .. .. ..$ total.variance: num 457
  .. ..$ umap:Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
  .. .. .. ..@ cell.embeddings           : num [1:2638, 1:2] -3.44 -7.94 -1.77 8.79 -1.51 ...
  .. .. .. .. ..- attr(*, "scaled:center")= num [1:2] 1.64 -1.3
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..$ : chr [1:2] "UMAP_1" "UMAP_2"
  .. .. .. ..@ feature.loadings          : num[0 , 0 ]
  .. .. .. ..@ feature.loadings.projected: num[0 , 0 ]
  .. .. .. ..@ assay.used                : chr "RNA"
  .. .. .. ..@ global                    : logi TRUE
  .. .. .. ..@ stdev                     : num(0)
  .. .. .. ..@ key                       : chr "UMAP_"
  .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
  .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ]
  .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ]
  .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ]
  .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ]
  .. .. .. ..@ misc                      : list()
  ..@ images      : list()
  ..@ project.name: chr "SeuratProject"
  ..@ misc        : list()
  ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
  .. ..$ : int [1:3] 4 1 3
  ..@ commands    :List of 7
  .. ..$ NormalizeData.RNA       :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "NormalizeData.RNA"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:21"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "NormalizeData(seurat)"
  .. .. .. ..@ params     :List of 5
  .. .. .. .. ..$ assay               : chr "RNA"
  .. .. .. .. ..$ normalization.method: chr "LogNormalize"
  .. .. .. .. ..$ scale.factor        : num 10000
  .. .. .. .. ..$ margin              : num 1
  .. .. .. .. ..$ verbose             : logi TRUE
  .. ..$ FindVariableFeatures.RNA:Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "FindVariableFeatures.RNA"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:22"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "FindVariableFeatures(seurat, nfeatures = 500)"
  .. .. .. ..@ params     :List of 12
  .. .. .. .. ..$ assay              : chr "RNA"
  .. .. .. .. ..$ selection.method   : chr "vst"
  .. .. .. .. ..$ loess.span         : num 0.3
  .. .. .. .. ..$ clip.max           : chr "auto"
  .. .. .. .. ..$ mean.function      :function (mat, display_progress)
  .. .. .. .. ..$ dispersion.function:function (mat, display_progress)
  .. .. .. .. ..$ num.bin            : num 20
  .. .. .. .. ..$ binning.method     : chr "equal_width"
  .. .. .. .. ..$ nfeatures          : num 500
  .. .. .. .. ..$ mean.cutoff        : num [1:2] 0.1 8
  .. .. .. .. ..$ dispersion.cutoff  : num [1:2] 1 Inf
  .. .. .. .. ..$ verbose            : logi TRUE
  .. ..$ ScaleData.RNA           :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "ScaleData.RNA"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:23"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr [1:2] "ScaleData(seurat, vars.to.regress = c(\"nFeature_RNA\", " "    \"percent.mt\"))"
  .. .. .. ..@ params     :List of 11
  .. .. .. .. ..$ features          : chr [1:500] "PPBP" "LYZ" "S100A9" "IGLL5" ...
  .. .. .. .. ..$ assay             : chr "RNA"
  .. .. .. .. ..$ vars.to.regress   : chr [1:2] "nFeature_RNA" "percent.mt"
  .. .. .. .. ..$ model.use         : chr "linear"
  .. .. .. .. ..$ use.umi           : logi FALSE
  .. .. .. .. ..$ do.scale          : logi TRUE
  .. .. .. .. ..$ do.center         : logi TRUE
  .. .. .. .. ..$ scale.max         : num 10
  .. .. .. .. ..$ block.size        : num 1000
  .. .. .. .. ..$ min.cells.to.block: num 2638
  .. .. .. .. ..$ verbose           : logi TRUE
  .. ..$ RunPCA.RNA              :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "RunPCA.RNA"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:25"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "RunPCA(seurat, npcs = 50)"
  .. .. .. ..@ params     :List of 10
  .. .. .. .. ..$ assay          : chr "RNA"
  .. .. .. .. ..$ npcs           : num 50
  .. .. .. .. ..$ rev.pca        : logi FALSE
  .. .. .. .. ..$ weight.by.var  : logi TRUE
  .. .. .. .. ..$ verbose        : logi TRUE
  .. .. .. .. ..$ ndims.print    : int [1:5] 1 2 3 4 5
  .. .. .. .. ..$ nfeatures.print: num 30
  .. .. .. .. ..$ reduction.name : chr "pca"
  .. .. .. .. ..$ reduction.key  : chr "PC_"
  .. .. .. .. ..$ seed.use       : num 42
  .. ..$ RunUMAP.RNA.pca         :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "RunUMAP.RNA.pca"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:35"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "RunUMAP(seurat, dims = 1:20)"
  .. .. .. ..@ params     :List of 26
  .. .. .. .. ..$ dims                : int [1:20] 1 2 3 4 5 6 7 8 9 10 ...
  .. .. .. .. ..$ reduction           : chr "pca"
  .. .. .. .. ..$ assay               : chr "RNA"
  .. .. .. .. ..$ slot                : chr "data"
  .. .. .. .. ..$ umap.method         : chr "uwot"
  .. .. .. .. ..$ return.model        : logi FALSE
  .. .. .. .. ..$ n.neighbors         : int 30
  .. .. .. .. ..$ n.components        : int 2
  .. .. .. .. ..$ metric              : chr "cosine"
  .. .. .. .. ..$ learning.rate       : num 1
  .. .. .. .. ..$ min.dist            : num 0.3
  .. .. .. .. ..$ spread              : num 1
  .. .. .. .. ..$ set.op.mix.ratio    : num 1
  .. .. .. .. ..$ local.connectivity  : int 1
  .. .. .. .. ..$ repulsion.strength  : num 1
  .. .. .. .. ..$ negative.sample.rate: int 5
  .. .. .. .. ..$ uwot.sgd            : logi FALSE
  .. .. .. .. ..$ seed.use            : int 42
  .. .. .. .. ..$ angular.rp.forest   : logi FALSE
  .. .. .. .. ..$ densmap             : logi FALSE
  .. .. .. .. ..$ dens.lambda         : num 2
  .. .. .. .. ..$ dens.frac           : num 0.3
  .. .. .. .. ..$ dens.var.shift      : num 0.1
  .. .. .. .. ..$ verbose             : logi TRUE
  .. .. .. .. ..$ reduction.name      : chr "umap"
  .. .. .. .. ..$ reduction.key       : chr "UMAP_"
  .. ..$ FindNeighbors.RNA.pca   :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "FindNeighbors.RNA.pca"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:35"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "FindNeighbors(seurat, dims = 1:20)"
  .. .. .. ..@ params     :List of 17
  .. .. .. .. ..$ reduction      : chr "pca"
  .. .. .. .. ..$ dims           : int [1:20] 1 2 3 4 5 6 7 8 9 10 ...
  .. .. .. .. ..$ assay          : chr "RNA"
  .. .. .. .. ..$ k.param        : num 20
  .. .. .. .. ..$ return.neighbor: logi FALSE
  .. .. .. .. ..$ compute.SNN    : logi TRUE
  .. .. .. .. ..$ prune.SNN      : num 0.0667
  .. .. .. .. ..$ nn.method      : chr "annoy"
  .. .. .. .. ..$ n.trees        : num 50
  .. .. .. .. ..$ annoy.metric   : chr "euclidean"
  .. .. .. .. ..$ nn.eps         : num 0
  .. .. .. .. ..$ verbose        : logi TRUE
  .. .. .. .. ..$ force.recalc   : logi FALSE
  .. .. .. .. ..$ do.plot        : logi FALSE
  .. .. .. .. ..$ graph.name     : chr [1:2] "RNA_nn" "RNA_snn"
  .. .. .. .. ..$ l2.norm        : logi FALSE
  .. .. .. .. ..$ cache.index    : logi FALSE
  .. ..$ FindClusters            :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "FindClusters"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:36"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "FindClusters(seurat, resolution = 1)"
  .. .. .. ..@ params     :List of 10
  .. .. .. .. ..$ graph.name      : chr "RNA_snn"
  .. .. .. .. ..$ modularity.fxn  : num 1
  .. .. .. .. ..$ resolution      : num 1
  .. .. .. .. ..$ method          : chr "matrix"
  .. .. .. .. ..$ algorithm       : num 1
  .. .. .. .. ..$ n.start         : num 10
  .. .. .. .. ..$ n.iter          : num 10
  .. .. .. .. ..$ random.seed     : num 0
  .. .. .. .. ..$ group.singletons: logi TRUE
  .. .. .. .. ..$ verbose         : logi TRUE
  ..@ tools       : list()

My objective is

Generate feature plot
Then label the cells based on only variable genes,right now I see all the genes I would like to only label based on top genes which I'm storing at top_gene but Im not sure how to incorporate it instead of every genes which I now see in my drop down menu
The after labelling cells i want to export the labelled cluster as csv file

What Now Im getting is html file which not what I want
Any help or suggestion is really helpful
UPDATE
The right function I have to use is this FetchData() but Im not sure how to pass the gene which i choose to label the cells and use the same gene as input for exporting cluster


text:
Protein fold topology search web-server
I'm finding a web server (or standalone program) that allows the upload of a custom PDB for protein fold topology search.
Do you know any of them?
More details:
Actually I started off with Foldseek but have a question in mind:
how do I know if I have a novel fold? (because similarity search will often return partially similar hits).
Then, I tried to think of topology, as an illustration:

which simplifies structures to spatial arrangements of secondary structures (or topology).
I want to know if there's any database of protein fold topologies and tools for topology similarity search.
ADD COMMENT • link • edit •

Have you tried foldseek?

I started off with foldseek but I need lower resolution search.

There are a few, but the main candidates for this would be :

CATH - http://www.cathdb.info/
SCOP - https://scop.mrc-lmb.cam.ac.uk/
Ptgl - http://ptgl.uni-frankfurt.de/

CATH and Ptgl have search services on their website, but not sure about SCOP.


text:
Comparing sequence similarity of two sets of peptide sequences?
What tools can we use to compare two sets of peptide sequences?

The biological question is I have some kmer sequences of some proteins and would like to know how similar are these kmer sequences compared with known MHC classs I and II peptides.
Therefore I need to compare two sets of peptides sequences.
The kmer and peptide sequences lengths range from 5-80aa.

Alternatively, I am also considering to use the full length of the proteins for comparison.

What kind of comparison? All against all? Pairs? Triplets? Are you expecting them all to be similar? Are you trying to group them into sets based on similarity? Blastp doesn't seem relevant at all here (and 11 is the default word size of blastn not blastp, but that's neither here nor there). Please [edit] your question and tell us what you are trying to do, what is the biological question you are trying to answer. We need to understand that in order to suggest approaches.

Hi @terdon, thanks for the suggestion and now edited the question.

Just a reminder to consider marking one of the questions as "answered". This is definitely helpful for all concerned.

Not sure what is the shortest peptide one can fish out with LAST out of small protein database (well, it starts with protein FASTA file), but you may check out LAST cookbook. This can use different PAM matrices so you have the proper similar amino acid matching. If that does not suit your needs there can be some more desperate route using string fuzzy matching i.e. with thefuzz python library and checking the best Levenshtein distance scores.

I would suggest focusing on short k-mers if you are using very small sequences. Simply quantify the similarity of the k-mer profile of each relevant pair of peptides. There are many tools for analyzing k-mer profiles, though many are targeted towards DNA so you should keep an eye out for that. Here is one python tool.
It is also not difficult to write a k-mer counting script if you decide that third party tools are not suitable.
Granted, you will need to evaluate that k-mer profile similarity has some relevance to your intuitive notion of "similarity", which it may not. It would clarify somewhat to know more about your desired application.

After more search, I go for blastp-short with MHC peptides as query and my full-length proteins as database. Please share if you think it has problems.


text:
How to write FASTA records using "Bio.SeqIO.write()"
Update:
Biopython document says that "Bio.SeqIO.FastaIO.FastaWriter" class is obsolete. Now my question becomes how to I use Bio.SeqIO.write() fucntion do achive the same thing I asked.
I use the following code to extract specific FASTA records from a file based on a list of FASTA IDs given in another file. The script works just fine. However, I want to able to write the results using "Bio.SeqIO.FastaIO.FastaWriter" class into a file instead of Python's print or write function. How can I do that in my last if statement?
from Bio import SeqIO
#import pprint
import textwrap

my_id_file     = open('sample_files2/fasta_id_records.txt','r')
my_fasta_file  = open('sample_files2/fasta_file.fasta','r')

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'
#pprint.pprint(my_dictionary)

for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:
        sequence = str(seq_record.seq)
        fasta_record = textwrap.fill(sequence, width=60)
        print(f">{seq_record.id}\n{fasta_record}")

my_fasta_file.close()
my_id_file.close()

I was considering modifying the title. I am glad you edited it @gringer. Thank you.

I figured it out:
from Bio import SeqIO
#import pprint
import textwrap

my_id_file     = open('sample_files2/fasta_id_records.txt','r')
my_fasta_file  = open('sample_files2/fasta_file.fasta','r')
result_file    = open('sample_files2/results.fasta', "w")

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'

#pprint.pprint(my_dictionary)
for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:
        #use this line to write standard fasta record using SeqIO.write
        SeqIO.write(seq_record, result_file, "fasta")

my_fasta_file.close()
my_id_file.close()
result_file.close()

I'm not 100% sure you can because that will be taken care within SeqIO.write(). The obvious seq_record.seq = textwrap.fill(str(seq_record.seq) , width=60) is unlikely to work.

The only way to do it is within a straight write to disk in you first code ... but (obviously) write rather than print.

it wraps to 60 by default in Biopython 1.8 at least, if saves as 'fasta-2line' uses a single line for the seq

my attempt at changing the way Biopython formats fasta file sequence width,
read into code to figure out standard behaviour:
from Bio import SeqIO, __version__
#import pprint
# import textwrap

print('Biopython version : ' , __version__)

# my_id_file     = open('sample_files2/fasta_id_records.txt','r')
# my_fasta_file  = open('sample_files2/fasta_file.fasta','r')
#result_file    = open('sample_files2/results.fasta', "w")
### CHANGED TO:

my_id_file     = open('fasta_id_records.txt','r')
my_fasta_file  = open('fasta_file.fasta','r')
result_file    = open('results.fasta', "w")

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'

#expected lenght of sequences in my fasta
length = 10

#pprint.pprint(my_dictionary)
for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:

        ########################################################################################
        #use this line to write standard fasta record using SeqIO.write

        # SeqIO.write(seq_record, result_file, "fasta-2line") # writes fasta in 2 line format

        # SeqIO.write(seq_record, result_file, "fasta") # writes fasta with width of 60 as standard biopython format

        ### NOT POSSIBLE !!!!!!!!!!!!
        # seq_record.seq = textwrap.fill(seq_record.seq , width=10) ### --> AttributeError: 'Seq' object has no attribute 'expandtabs'

        # print(seq_record.seq, type(seq_record.seq), dir(seq_record.seq))

        sequence = str(seq_record.seq)

        result_file.write('>'+seq_record.id+'\n')

        while len(sequence) > 0:
            result_file.write(sequence[:length]+'\n')
            sequence = sequence[length:]
        result_file.write('\n')

        ### if you try to add EOL inside seq_record.seq of course you get
        ### TypeError: SeqRecord (id=seq1) has an invalid sequence.
        # SeqIO.write(seq_record, result_file, "fasta") # writes fasta with sequence width as 60   !!!! TypeError: SeqRecord (id=seq1) has an invalid sequence.

my_fasta_file.close()
my_id_file.close()
result_file.close()

input files :
fasta_id_records.txt:
seq1
seq2
seq3

 fasta_file.fasta:
>seq0
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>seq1
TTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTT
>seq3
CCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGGGGGGGG
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

output file:
v :
>seq1
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTT

>seq3
CCCCCCCCCC
CCCCCCCGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGCCCCC
CCCCCCCCCC
CCCCCCCCCC
CCCCCCCCCC
CCCCCCCCCG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
G

You could instead pass in a generator that yields SeqRecords. This avoids storing the entire list of  sequences in memory and is usually the way to go assuming that you need to handle large files containing thousands of records. This example use a generator expression, which is a bit like a list comprehension, but instead of creating a list, it creates a generator:
import argparse
import contextlib
import pathlib
import sys

from Bio import SeqIO

def get_argument_parser():

    parser = argparse.ArgumentParser()

    parser.add_argument(
        "fasta_file",
        type=pathlib.Path,
        metavar='FILE',
        help="Extract sequence records from FILE",
    )
    parser.add_argument(
        "ids_file",
        type=pathlib.Path,
        metavar='FILE',
        help="Sequence record IDs/names to extract, one per line",
    )
    parser.add_argument(
        "-o",
        "--output",
        type=pathlib.Path,
        metavar='FILE',
        default="-",
        help="Write extracted sequences to FILE (default: stdout)",
    )
    parser.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )

    return parser

def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:

        ids_file = stack.enter_context(open(args.ids_file))
        fasta_file = stack.enter_context(open(args.fasta_file))

        if args.output != '-':
            out = stack.enter_context(open(args.output, 'w' if args.force else 'x'))
        else:
            out = sys.stdout

        record_ids = {line.rstrip('\n') for line in ids_file}

        fasta = SeqIO.parse(fasta_file, "fasta")
        records = (rec for rec in fasta if rec.id in record_ids)

        SeqIO.write(records, out, "fasta")

if __name__ == "__main__":
    main()

cool never eard about ExitStack till now https://docs.python.org/3/library/contextlib.html#contextlib.ExitStack

A better description about what it does is here https://www.rath.org/on-the-beauty-of-pythons-exitstack.html @pippo1980

I think I need to ask another question on this code. It's certainly professional code. What I'm wondering is how the if-condition tuple comprehension is a generator because thats the "workhorse" of the code.

@M__ The parentheses just make it a generator expression. Brackets would have made it a list comprehension: https://peps.python.org/pep-0289/

@Steve, gave me a code to study.


text:
How to visualize Similarity score matrix with CoNet app of Cytoscape
i got a similarity score matrix by using CCREPE package on R. I would like to visualize the co-occurence by using the app CoNet from Cytoscape. Everytime i upload the .TSV matrix and .TSV p-values, i got an error.
This is the similarity matrix:

This is the p-value matrix:

This is how i upload the tho files in CoNet:

This is the error i get once i tap "go":

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Please show some of your input data. The first few lines of the .TSV matrix and .TSV p-values would be helpful.

thanks for your reply. I posted more pictures

My guess is that CoNet doesn't like NaN entries. Try changing those to 0 and see if it works.


text:
How can I import xgmml cytoscape networks into a new session file such that style formatting can be applied?
I am attempting to consolidate several networks from across several Cytoscape sessions into a new session file for publication purposes. The goal is to export my list of individual networks from the respective source *.cys sessions and import them as fully editable networks with layout and style retained. I exported the original networks as *.xgmml files (with style info) and imported them into a new *.cys session. Style changes (e.g. border color, etc) can be modified for the first imported network, but all subsequent imported networks appear to be locked to style changes. (This happens irrespective of import to the same network collection or a new network collection.) Is there a fix / workaround?

Yes, as you noted, when you bring in an XGMML file, the style properties often come in a "locked".  Actually, the first one should also have been "locked".  You can remove the locked attributes by selecting the node or all nodes and then right-clicking on the bypass column and doing a "remove bypass".  My suggestion, however, would be to export the styles "File->Export->Styles to File" and just name them different things.  I think that will get you more of what you want.


text:
Extracting metadata for a list of BOLD BINs
When I view the Public Data Portal page for a specific BOLD BIN, I get some useful metadata, including BIN Details (the image below is from this BIN):

I would like to be able to compile and download this data for a list of BINs, like for example, all the BINs associated with the Family Braconidae. I could not find a way to do this through the BOLD website. Is there a tool or an API plugin that would allow me to do this? This is not strictly a bioinformatics question, so if there is a stackexchange that is better suited I can move it there.
Edit:
To be specific, I would like to know the minimum, maximum, and mean distance within a given set of BINs, and the distance between those BINs their nearest neighbours, as is displayed in the image above.
Thanks,
Y.

What do you need exactly? Yes, the site has APIs, see https://boldsystems.org/index.php/api_home but you need to show us exactly what fields you want in order for us to be able to answer.

Thanks @terdon , I have edited the question to include specific details.


text:
Cannot invoke "Object.toString()" but Uniprot IDs are there
I'm trying to stringify my gene network in Cytoscape and have designated a column with uniprot IDs for the basis of the search. Not every gene has a uniprot ID but I have selected to retain genes without them. For some reason I am unable to execute this function and receive  the feedback:
"Stringify network
 Cannot invoke "Object.toString()" because the return value of "java.util.List.get(int)" is null"
I have verified the uniprot IDs are correct and exist with singular matches to the database, could anyone give insight into how I might overcome this issue?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Could you provide a backtrace when you get this error?  You can get it by running Cytoscape from a command shell, but you can also get provide the backtrace that shows up in your CytoscapeConfiguration/3/framework-cytoscape.log  That will help us figure out exactly where the problem is.


text:
Differential gene expression (aiming to create heatmap) for raw RNAseq count data from different publications - how to?
I'm new to bioinformatics and RNAseq analysis, I'll try my best to explain the question!
I have raw counts RNAseq data from 3 different publications (so 3 different datasets) that all had similar methods to produce their data. I will be comparing the gene exp of Cell_type_A vs Cell_type_B in each of the datasets, so I will be doing differential gene expression on these 2 cell types. The 2 cell types will be my "conditions", if that makes sense. To do this, I used R to apply DESeq2 on each dataset (independently) to get the log2fold change of each gene for the differential expression of Cell_type_A vs Cell_type_B. I filtered for genes with p-adj < 0.05. Now that I have significant log2fold changes for each dataset, how can I compare these results in a heatmap? The type of heatmap I'm going for is 3 columns: one for each dataset. The heatmap will just display the log2foldchange values, like how a volcano plot would. Positive values (upregulated) = red, negative values (downregulated) = blue.
I was thinking of making a new dataset with 4 columns: gene, log2fs of dataset 1, log2fs of dataset 2, log2fs of dataset 3. And from there I can make a heatmap out of the new dataset: so my heatmap will have 3 columns of data and the genes as the row labels. But that seems like a cheap solution to me. Is it a valid way of going through with this? I really appreciate all kinds of help, I've been stuck on this problem for quite a while!

The question lacks any background information, e.g. what experiment is being performed and what coding can you do, presumably R maybe Python? You should perform a Volcano plot rather than a heat map. If I get time I will forward an answer.

Thank you for these suggestions, I will edit my post!

There's a lot to unpack here; it should really be in a few separate questions. I'm just going to concentrate on the data preparation side of things....
Given that you're using DESeq2, it would be better to prepare a single count table with all your data (genes / features in rows, samples in columns), and add any complicating factors into the statistical model.
Given that you're working with multiple datasets, one of those complicating factors would be the dataset that the data came from. So your meta data frame (e.g. meta.df) should look something like this:
Dataset Condition SampleID
Pub1    Cell_A    NCR00102
Pub1    Cell_A    NCR00103
Pub1    Cell_A    NCR00104
Pub1    Cell_B    NCR00112
Pub1    Cell_B    NCR00113
Pub1    Cell_B    NCR00114
Pub1    Cell_B    NCR00115
Pub1    Cell_B    NCR00118
Pub1    Cell_B    NCR00120
Pub2    Cell_A    NBS00233
Pub2    Cell_A    NBS00234
Pub2    Cell_A    NBS00235
Pub2    Cell_A    NBS00236
Pub2    Cell_A    NBS00238
Pub2    Cell_B    NBS00030
Pub2    Cell_B    NBS00031
Pub2    Cell_B    NBS00037
Pub2    Cell_B    NBS00042
Pub3    Cell_A    HGZ00780
Pub3    Cell_A    HGZ00878
Pub3    Cell_A    HGZ01023
Pub3    Cell_A    HGZ02053
Pub3    Cell_A    HGZ03950
Pub3    Cell_A    HGZ04059
Pub3    Cell_A    HGZ05069
Pub3    Cell_B    HGZ02042
Pub3    Cell_B    HGZ03056

Then the DESeq2 command would look something like this:
dds <- DESeqDataSetFromMatrix(countData = counts.mat,
                              colData = meta.df,
                              design= ~ Dataset + Condition)
res <- results(dds, name="Condition_Cell_B_vs_Cell_A")

Information about doing this can be found in the DESeq2 manual. I'd recommend reading at least the quick start section; many of the answers to my own questions about differential analysis can be found in that manual:
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#quick-start


text:
What parameter and threshold to choose to get significant pathways from MetaboAnlyst pathway analysis results?
MetaboAnalyst pathway analysis results give a total of 4 different p-values viz, raw p-value, -log10 pvalue, Holm-Bonferroni adjusted p-value and the FDR value. And an impact score.
Which value should be considered in my case (above image) to get significant pathways? And what should be the threshold cut off?

It is best recommended not to upload a picture of text as it becomes hard to read and reproduce. Instead, create a table using MarkDown (As shown here - https://meta.stackexchange.com/questions/356997/new-feature-table-support). If the image is from a journal/article, please list the article and table number so people can access it and help you better.

There are a number of factors to consider when doing a pathway Over Representation Analysis (ORA), in particular for metabolomics. I recommend an overview such as https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8448349/
This covers concepts such as correcting for multiple hypotheses (as well as other considerations). You should always do this to avoid false positives (i.e don’t use the raw p values). Note that whichever correction method you use for your results, none of the pathways are significant, so the specifics are somewhat moot for these results. But even when you do have significant corrected p-values, overrepresented should be treated as hypotheses.


text:
Deseq2 model formula and longitudinal experimatal (or time-series) designs
At hand I have RNA expression data of 6 organs from mice with treated and untreated diets, collected at consecutive time points, and I wonder how to design an analysis accordingly in DESeq2. I am rather new to applying DESeq2.
Dietary treatment was started at 0 minutes, and now RNA sequencing data data is available for each of the time points at 0, 30, 60, 120, 240, and 480 minutes, and for each dietary treatment (TREAT) and control (CNTRL) (So, in total I have about 72 quadruplicates from 6 organs x 6 time points x 2 diets (TREAT/CNTRL)) .
I was wondering what the correct way is to obtain information about the differentielly  expressed genes (DEGS) between:

treatment and control - at each time point
between each consecutive time point - within each diet (treatment and control)

The two possible solution I am consided are:

As I encoded already as per Fig 1 below: I divided my data (for each organ) into sub-sets (CNTRL vs. TREAT at each time point, and TREAT or CNTL for each consecutive time point) and calculated DEGs individually - resulting in 16 distinct DEG lists (Fig 1 below). The DESeq2 model formulae had only one term - time, or diet. To me this is a correct way - as the relationship between two time points or diets can adequately be described as a line - i.e. 1 ~ diet at each time-point can only be meaningfully described as a (regression) line, and 1 ~ time for 2 time points with each treatment can only be meaningfully described as a line. But then I assume I can't compare effect sizes.

I am wondering whether it would be more appropriate if I wanted to compare results to include both time and ant diet into the model formulae, and then pull out 16 different results list by supplying 16 different contrast - defined following the schematic in Fig. 1. To me that would make the DEG lists comparable to each other, but it appears incorrect because the real observations across all time points could likely hardly be solved by linear regression - in other words the relationship 1 ~ diet * time across all genes is not a very large collection of lines with different slopes, but rather a collection of curves that can't be solved with linear regression.

I'd be happy to receive some feedback or thoughts how to approach this analysis In essence the one question I am entertaining is: Would you also choose approach 1, or rather 2, or a different approach entirely?
Kind regards, Paul

I appreciate the level of detail in the question, however please edit and focus the post on a single question. You can then always ask follow-up questions. The site works on a one question - one answer format.

Thank you for the hint, I have edited my post to only contain one question.

There are a few ways to approach your question, as you're partly asking for an opinion on the best method. In my view, to analyze changes in gene expression over time between treatment and control groups, it’s obviously crucial to focus on the interaction between these variables in your model.
I’d recommend setting up your DESeq2 model to include both time and treatment, along with their interaction term: ~ organ + diet + time + diet:time.
This setup allows you to focus on the interaction term (diet:time), which will highlight genes where the expression patterns over time differ between TREAT and CNTRL.
Once you’ve looked at the interaction term, you can still go ahead and examine the TREAT vs. CNTRL contrasts at each specific time point to see when those differences are most significant.
As for your concern about the linear regression not capturing non-linear relationships—it's a valid point. However, DESeq2 is generally "robust enough" for RNA-seq data. It really depends if you think that tradeoff is worth it because it allows you to compare estimates of variability and effect sizes across different conditions.

Thank you very much for your reply I appreciate the exchange.


text:
How to adjust QIIME DADA2 parameters to avoid eliminating real sequences
I am working with sequence reads (acquired through MiSeq platform by Illumina) from a viral inoculum. I have already used FastP to merge paired-end reads and cut 5' and 3' ends to remove adapters and low Phred scores (threshold = 32). The sample should exhibit a genetic diversity that I will label as X. X is defined as the "true theoretical" average sequence identity of all of the reads that should be present in the inoculum (after removing reads with PCR errors and hypermutation and without removing any of the "real" sequences).
The overall genetic diversity of my inoculum is too high for what I expect given the natural mutation rate of this particular virus (i.e. measured diversity >> X). Previous attempts at denoising or eliminating artefacts still has not lowered the genetic diversity to an acceptable level. Lately, I have been using the DADA2 plugin for QIIME2 (version qiime2-amplicon-2024.5) to denoise my reads, as this has produced the most success in eliminating artefacts. However, it looks like it's been working too well because it has lowered the overall genetic diversity to a level lower than X. This would suggest that DADA2, with my given params/args/options, is not only eliminating artefacts but some of the real sequences as well.
Is there a set of QIIME DADA2 params/args/options that would help avoid eliminating real sequences during denoising and would achieve a final variant diversity as close to X as possible? Here is the code that I'm using for denoising:
qiime dada2 denoise-single \
--i-demultiplexed-seqs FILE.qza \
--p-trunc-len 0 \
--p-max-ee 2 \
--p-pooling-method "independent" \
--p-chimera-method "none" \
--p-min-fold-parent-over-abundance 45 \
--p-allow-one-off \
--p-n-reads-learn 1000000 \
--p-no-retain-all-samples \
--o-table FILE-table.qza \
--o-representative-sequences FILE-seqs.qza \
--o-denoising-stats FILE-stats.qza

Thank you for any help you can give me!

I am not sure that DADA2 is the right tool for this, though it might help to hear more about your experimental design. I'd recommend reading up on it or looking at materials like this.
It seems that the mutation rate is relevant to your analysis, such that you are calibrating expected diversity $X$ based on it. Depending on the expected frequency distribution of your variants, DADA2 might be filtering out most of the variation that you're interested in as sequencing errors.
See this note on the ASV approach used in DADA2 (reference):

ASVs are inferred by a de novo process in which biological sequences are discriminated from errors on the basis of, in part, the expectation that biological sequences are more likely to be repeatedly observed than are error-containing sequences. As a result, ASV inference cannot be performed independently on each read—the smallest unit of data from which ASVs can be inferred is a sample. However, unlike de novo OTUs, ASVs are consistent labels because ASVs represent a biological reality that exists outside of the data being analyzed: the DNA sequence of the assayed organism. Thus, ASVs inferred independently from different studies or different samples can be validly compared.

In brief, you will only be able to infer the presence of mutants that occur at a rate detectably higher than the rate of sequencing error.
While I am not a DADA2 expert, what DADA2 is expecting is more along the lines of multiple modes in the sequence distribution- each of which is surrounded by a small scattering of sequencing errors resulting in variants nearly identical to the parent mode. It then uses some statistics to evaluate how likely each variant is to be an error off of the nearest mode. So variants that are very close to an existing mode are likely to be lumped into that mode, unless they are sufficiently high abundance to justify establishing a new "species" mode. This has led to a lot of success for DADA2 in accurately identifying variation in samples of relatively diverse samples with somewhat distantly related organisms each at some meaningfully high relative abundance (>=0.1%).
See for example this figure from the original DADA2 paper:

Every point is a variant identified by DADA2. you'll notice that the variant frequency has to be quite high (~0.1%-1%) to identify distinct variants from sequencing error at a Hamming distance of less than 10 or so from the nearest other variant. Unless you expect your mutant frequencies to be in that range, you are unlikely to detect them. Maybe they've gotten better over the years, but I expect that there are some fundamental statistical issues here for this approach that mean they can't get much better.
If you really want to use DADA2, any parameter that loosens the statistical thresholds will make it better, but I'd guess you've already tried that. I don't know the QIIME parameters, but if you look at dada2 docs, you will see that the dada() function has a number of error parameters that you could play with to try to tune down the probability of errors to increase the likelihood of recovering variants. See the learnErrors() function that you might be able to play with to get reasonable values for this, if you had e.g. a pure reference virus sequenced sample. Note that these functions are for denoising a diverse metagenomic community, not a low-diversity viral mutation experiment! Finally, see setDadaOpt() as ways of globally controlling various thresholds.
I would recommend doing this within the R package rather than through QIIME, to give you better control. If you can't do that, then maybe re-estimating your FILE-stats.qza or your p-min-fold-parent-over-abundance parameter might let you loosen thresholds a bit (though it looks like you are already doing this last somewhat).
outside forums
see also some similar questions on forums.


text:
Finding transcriptome matches for thousands of 21-mers
I am on a project to find sequence matches in mRNA databases across 20+ taxa. I was presented an EPA memorandum that outlines a method I was requested to replicate. I am not sure it is the best method and the memorandum didn't go into much detail. I played around a bit and could use a bit of advice.
Setup:
I have an unspecified number (I haven't been told yet) of dsRNA segments approximately 300bp in length. I need to match these to the human transcriptome as well as over 20 other taxa.
Criteria:
For each 300bp dsRNA, I am to find mRNA in the taxa that have 14 or more matches within a 21bp window. Then I sort the data by taxon, transcripts matched, and the annotation for the matched mRNA.
Approach (this is where I have questions):
The EPA memorandum says the Burrows-Wheeler Aligner (BWA) was used to align a 21-mer sliding window along the target transcriptomes to look for matches of 14 or greater within the window. The PI said to create all 21-mers using a sliding window along the dsRNA sequence. Easy enough.
Here are my questions:

Is BWA the best approach to use? I've never used BWA MEM for anything so small.
How should I set the parameters for the BWA for this case? The defaults are inadequate, but I'm just taking stabs in the dark to see what falls out. So far, I have adjusted:

Minimum seed length (-k) down to 3
band width (-w) down to 7
ignore alignment scores lower than (-T) range from 1 to 21
gap open penalty (-O) between 1 and 6
mismatch penalty (-B) between 1 and 4

Why do I see a Bitwise Flag of 0? In adjusting the parameters, the resulting SAM will contain matches where the Bitwise Flag is 0. This seems like nonsense to me, suggesting that I may be on the wrong track.

Sample Execution:
./bwa mem -k 5 -B 1 -O2 -T 5 ../ncbi_dataset/GCF000001405.40.rna.fna ../seqA.fasta | gzip -3 > ../bwa_results/aln_seqA.sam.gz
Bitwise Flag == 0?
seqA_331_352   16      NM_014873.3     6590    0       6S15M   *       0       0       GATCGGTGTAAATCCCATATC   *       NM:i:1  MD:Z:14A0       AS:i:14 XS:i:14
seqA_332_353   16      XM_011510229.4  7276    0       7S14M   *       0       0       TGATCGGTGTAAATCCCATAT   *       NM:i:0  MD:Z:14 AS:i:14 XS:i:14
seqA_333_354   0       XM_017008212.3  1680    0       7S14M   *       0       0       TATGGGATTTACACCGATCAA   *       NM:i:0  MD:Z:14 AS:i:14 XS:i:13
seqA_334_355   0       XM_017008212.3  1680    0       6S15M   *       0       0       ATGGGATTTACACCGATCAAC   *       NM:i:1  MD:Z:14A0       AS:i:14 XS:i:13
seqA_335_356   0       XM_017008212.3  1680    0       5S14M2S *       0       0       TGGGATTTACACCGATCAACT   *       NM:i:0  MD:Z:14 AS:i:14 XS:i:13```

What is it you consider small for bwa? If anything, 300bp is on the long side. And you only have single, individual sequences, right? Not multiple reads with the same sequence? What am I missing? And have you considered using blast instead? It's the kind of thing blast is made for.

The 300bp sequence isn't fed into BWA. Rather, 21-mers are generated using a sliding window. It is the 21-mers that are put into BWA.

Depending on the input 300bp query sequences (A, B, C and so on) you may get identical 21-kmers derived from the same (say A RNA) or several (A, B) sequences. To skip the repetitive searches one can cluster at 100% identity 21-kmers keeping track of the origins using i.e. vsearch
If you must use bwa mem try faster bwamem2.
But before plunging into using bwa or blast I suggest doing some test runs with vsearch or LAST to compare the results/do the benchmarks.


text:
Nextflow Best Practice When Collecting a Large Number of Files in One Step
Consider a hypothetical pipeline shown below, where many files need to be collected at one stage in the pipeline:

It seems like the idiomatic way of accomplishing this is to pass Process B a list of filenames, which then passes those file names to the tool in question:
workflow {
    channel_with_many_files = channel.of(...)

    PROCESS_A(channel_with_many_files)

    all_files_together = PROCESS_A.out.collect()

    PROCESS_B(all_files_together)
}

And then in the process B definition:
process PROCESS_B{
    input:
    path(all_files_list)

    output:
    path(single_output_file)

    script:

    """
    ./mytool ${all_files_list}
    """
}

The problem with this, in my opinion, is that when n (number of files) grows to be very large (say > 2000), the command that is generated in PROCESS_B is unweidly and a bit messy.
Is this the way to deal with such a pattern in Nextflow, or is there a better way?

A better way would be to stage the files in a sub-directory, which would allow you to simply pass in a glob pattern to your tool. This can be achieved using the path input type and a string containing the * or ? wildcards to control the names of staged files1. Here, the files are named as the source files in a sub-directory called 'staged':
process PROCESS_B{

    input:
    path 'staged/*'

    script:
    """
    mytool staged/*
    """
}

Note that you can also use the stageAs option for this, if a variable is somehow useful:
process PROCESS_B{

    input:
    path all_files_list, stageAs: 'staged/*'

    script:
    """
    echo "Staged ${all_files_list.size()} files..."
    mytool staged/*
    """
}

Oh man, I can't believe I missed this in the docs. This is very elegant. Thank you!


text:
KeggScape Network Import Issue
Whenever I try to import an xgml file in cytoscape the file shows gibberish nodes like some form of code for example eco 234 all the nodes have similar codes like this instead of the actual name of the molecules in my network. How do I fix this?
Update: Hello there as of 26/10/24, there has been an update for Cytoscape to the version 3.10.3. If you are updated and all including having the updated keggscape app, try uninstalling keggscape and reinstalling the latest version and I can now confirm that the pathways are showing up as expected thankyou for all the help.
ANOTHER UPDATE 10/30/2024: SO the networks are showing up as expected individually but when I merge them the same problem is occuring. Please help again.

Keggscape is not correctly installing the style.  We're aware of the problem and are working on it...and have now just fixed it.
Please install the latest version.

Can you kindly described the versions where the bug was discovered and the version where it was fixed? The reason is that each post exists long after the posting date and thus will be a reference to others.


text:
What are some good protein databases to use? Trying to find neuraminidase sequences for influenza
I am going to study neuraminidase and hemaglutinin, proteins of interest for influenza virus vaccines. My first step is going to be to find all the protein sequences I can, make databases of my own in order to study them. I started using uniprot mostly and looking for similar proteins with BLAST. My main question is actually a string of related questions that I am going to leave here:
Does Uniprot have all the proteins that are known? If not, is there such a database? If not, do I combine results from multiple databases? What would other good databases like Uniprot be?
I have other questions like: How do i make sure i dont leave out some proteins? Is it worth thinking about this or is it inevitable? etc. but I dont think asking everything here would be efficient.
Thank you for your time! (And sorry for any dumb question)
EDIT(context): I am looking for neuraminidase sequences from influenza, the final purpose is to find epitopes and the corresponding antibody structures to interact with those epitopes. As for the specific analyses that I will do, I am in a very early stage and I am not sure of that myself as I am still learning. As for blast, I am getting acquainted with it now, i used it to find similar sequences to the ones I selected, though I am also unsure about what to select and how to properly use it. I am reading the documentation and trying to figure those things out at the moment but this is the general context.

Yes, UniProt is the most comprehensive database of protein sequences we have. Or, at any rate, the most comprehensive one I know. However, if you're looking for neuraminidases, don't try and find them using blast. While blast can indeed be part of a pipeline designed to find homologous sequences, it isn't enough by itself and just running blast to get similar sequences is not a good way of finding actual homologs if you don't know what you're doing, apply correct filters etc.
Instead, if you have your species of interest, you can look up the proteins in UniProt or simply NCBI's RefSeq databases. You also have various databases that host collections of homologous sequences. Don't reinvent the wheel. You might want to post a new question, explaining exactly what kind of sequences you need and asking for help on obtaining them. But really, don't try and just blast your way to a set of homologs, that is quite a complicated job to get right.


text:
When calculating Faith pd with picante package in R species richness appears very low for most samples but that is not true
I am trying to compute Faith's pd for the first time using the picante package in R.
My code is the following:
library(picante)
library(tidyverse)

set.seed(19990409)

setwd("~/Desktop/Urbanization/Bacteria_QIIME2/Urb_Grad")

# loading a tsv list that has all OTUs as the first column, then all samples with the subsequent abundance of the OTU and taxonomy
#  data has been rarefied to match the sample with the fewest observed features

tab_samp_feat <- read_tsv("feature_and_taxonomy.txt",col_types = list(OTU_ID=col_character(),
                                                                     taxonomy=col_character(),
                                                                     .default=col_number()))

tab_samp_feat <- tab_samp_feat %>% select(-taxonomy) %>%
  rotate_df()
colnames(tab_samp_feat) <- tab_samp_feat[1,]
tab_samp_feat <- tab_samp_feat[-1,]

metadata <- read_tsv("Metadata.csv", col_types = cols(Samples = col_character(),
                                                                                Sampling_season = col_character(),
                                                                                Site = col_character(),
                                                                                .default = col_number())) %>%
  select(c("Samples", "Sampling_season", "Site")) %>% as.matrix()
rownames(metadata) <- metadata[,1]
metadata <- metadata[,-1]

# rooted tree exported from QIIME2
phy <- read.tree("Exported/tree.nwk")

combined <- match.phylo.comm(phy, tab_samp_feat)

phy <- combined$phy
dat <- combined$comm
metadata <- metadata[match(rownames(dat), rownames(metadata)), ]
all.equal(rownames(dat), rownames(metadata))

dat.pd <- pd(dat, phy)
head(dat.pd)
boxplot(dat.pd$PD ~as.data.frame(metadata)$Site, xlab = "Site", ylab = "PD")

Everything seems to work alright, in the sense of all OTUs in my samples existing in the tree and the order of the metadata samples being the same as the order of the data. But the output is very strange:
PD results:
The PD metric is very low, but I think that is a consequence of the very small SR number for most samples. However I believe that is wrong as for example, even when looking at the first sample M-N1 there are way more than two OTUs present.
Screenshot of more than 2 features being present in the M-N1 sample
I was expecting the SR column to reflect the real presence/absence data of the features in my sample and for SR and PD to have more balanced values.
I would be so grateful if anyone could help me realize what it is that I am doing wrong.
Thank you!
Update: I tried using QIIME2 to visualize the Faith pd metric and the results look like the following:

The question has been given a tag of "microbiology". What species groups are being analysed and any biogeographic or sample information. Is this metagenomics data?
Basically, the question cannot be answered without surrounding biological information, because its not clear whether "high" or "medium" phylogenetic diversity should be expected.

I assume the OP has correctly rooted all trees in the analysis consistently. If this has not happened then that could be problematic and explain the discrepancy.
Otherwise ....

Basically I suspect this data is single locus metagenomics. Build a tree and quantitate the phylogenetic distance as a metric.
The pd scores between the figure and the individual samples are not comparable, because the figure appears to represent grouped samples. The values from individual samples need some definition about what they represent. More generally an additional layer of analysis is required like Kraken2 to underpin what the diversity is and whether it is representative of the sample.
Basically, if this is all the information provided, the pd scores are not necessarily low, its just comparing grouped data to individual samples. Grouped data will inevitably cause a large shift in genetic diversity represented by pd.
Thus two samples of pd = 3 could generate a pd of 40 when combined, but it depends on the species composition of the samples.


text:
Requesting help with Edge Weighted Graph Creation (Cytoscape)
I've used WGCNA to calculate adjacency and TOM for a set of 20k genes and particularly extracted adjacency and TOM for one gene relative to the 20k other genes. I would like two create a graph on cytoscape with the one gene as the central gene and the other genes to flank around it with respective to their adjacency or TOM with the gene.
Using the cytoscape app and importing the edges + nodes data (edge data loaded as network first and then node data loaded as a table into he network), I've been able to do this but I keep encountering an issue with accuracy of the mapping. I've used the edge weighted spring embedded layout but I get visualizations with some inaccuracy (i.e some genes with a higher adjacency may further away from the central gene than genes with lower adjacency even though the rest of the graph follows the intended trend).
There is a problem with randomness, when choosing the edge weighted spring embedded layout on the basis of adjacency as edge weight, the layout may randomly arrange genes with higher adjacency closer to the gene sometimes and vice versa other times. When I try to manually change settings in layout settings and apply the layout based on those settings, the applied layout won't even reflect the selected column for "edge column that contains weights", the distribution of genes around the central gene is random. I have colored nodes on the basis of pearson correlation and their shape height and width to be reflective of their average expression - if that has any impact on layout. All the weights (adjacency, TOM, etc.) have been loaded in as floats.
Would anyone happen to have any input regarding the randomness I'm seeing, it's unsettling that there are issues with accuracy and a lack of irreproducibility. I am using the latest cytoscape version (3.10.3). Thank you!


text:
BCFTOOLS: command does not work correctly. Can you help me?
Cross-posted on biostars

So, I used the command bcftools view -v indels SRZ189891_722g.990.SNP.INDEL.chrAll.vcf.gz -Oz -o indels.vcf.gzto create a file with only the indels. When I checked if the number of indels in the file SRZ189891_722g.990.SNP.INDEL.chrAll.vcf.gz was the same as in indels.vcf.gz, the results matched. However, I wanted to verify if any SNPs had been filtered out as well, so I used the command bcftools view -v snps indels.vcf.gz|grep -v -c '^#'. The result was that there were SNPs in the file indels.vcf.gz, when there should only be indels.
What could have happened? What can I do to have only the indels?
Edit: here are my data
chrX    77397836        .       CAG     CAGAGAG,CAGAG,GAG,C,CAGAGAGAGAG,CAGAGAGAGAGAGAG
chrX    77398227        .       T       A,TA
chrX    77398233        .       TTTATTTAA       ATTATTTAA,T
chrX    77398237        .       TTTAATTAA       TTTAA,ATTAATTAA,T
chrX    77398241        .       A       T,AT,ATTTATTTAT
chrX    77400094        .       C       T,CA
chrX    77401622        .       CA      C,AA
chrX    77404139        .       CG      C,TG
chrX    77406406        .       TCC     TCCC,CCC,TCCCC,T
chrX    77411156        .       CAAATAAATAAAT   CAAAT,CAAATAAAT,CAAATAAATAAATAAAT,TAAATAAATAAAT,C

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

The first step of any analysis should be to actually check your data. I tried your commands on a vcf file I had, and the answer was obvious:
$ bcftools view -v indels input.vcf.gz -Oz -o indels.vcf.gz
$ bcftools view -v snps indels.vcf.gz | tail -n3 | cut -f 1-5
22  49591397    .   T   A,TGTAG
22  49839228    .   G   A,GAGAA
22  50602312    .   AT  TT,A

As you can see, in my case and most likely in yours, what happens is that you have one line describing multiple variants. This is relatively common. You are looking at cases where the snp is heterozygous and the sample has an indel at the same position as the snp. Or you might be dealing with cancer data or mosaicism or non diploid organisms, there are loads of different reasons this might happen.
If you want only snps even in cases where you have both a snp and an indel in the same position, one approach is to first split all multi-allelic sites, so that two variants in the same position are now shown on separate lines, and then select snps only.
So first use bcftools norm -m -any to split. The -m option is:
    -m, --multiallelics -|+TYPE     Split multiallelics (-) or join biallelics (+), type: snps|indels|both|any [both]

So -m -any will split all multiallelic sites. Now you can use view to get indels only:
bcftools norm -m -any input.vcf.gz -Oz |
  bcftools view -v indels -Oz > indels.vcf.gz

If you now check the new indels.vcf.gz file, it should contain only indels:
$ bcftools view -v snps indels.vcf.gz|grep -v -c '^#'
0

And looking at one of the examples from before, we can see that only the indel has been kept:
$ zgrep 49591397 indels.vcf.gz | cut -f 1-5
22  49591397    .   T   TGTAG

While the original had both:
$ zgrep 49591397 input.vcf.gz | cut -f 1-5
22  49591397    .   T   A,TGTAG

Sorry for the delay in responding. I edited my comment and added the data. I used the command bcftools view -v snps indels.vcf.gz | tail -n10 | cut -f 1-5. Hope this helps. I'm new to bioinformatics and know almost nothing about it.

@JhonatanHenriqueRossiniLcio see updated answer for how to keep only indels in the end.

I have a question: do I first use the command bcftools norm -m -any input.vcf.gz -Oz and then the command bcftools view -v indels -Oz > indels.vcf.gz or do I use the command bcftools norm -m -any input.vcf.gz -Oz | bcftools view -v indels -Oz > indels.vcf.gz all at once? Thanks for helping me

@JhonatanHenriqueRossiniLcio you run them both together. That's why I have them both together in the answer.

Thanks for helping me out! So, I used the command bcftools norm -m -any input.vcf.gz -Oz | bcftools view -v indels -Oz > indels.vcf.gz and it gave me a file. When I used the command bcftools view -v snps indels.vcf.gz|grep -v -c '^#' to see how many snps there were, there were none. However, when I used the command bcftools view -v indels indelss.vcf.gz|grep -v -c '^#' to see if it had filtered out the indels correctly, it gave me a value of 26067197, when there were originally 18722319. What could have happened?

@JhonatanHenriqueRossiniLcio I don't know. Have you tried looking? We don't have access to your data, you will need to put some effort yourself. Find specific indels present in one file and absent in the others and see if you can figure it out. If not, please post a new question since this site is only for Q&A and not discussion.

okay. thank you very much for your help


text:
Metabolite function from DNA sequences
I have an aligned metagenome fasta file. I want to see the relation between the aligned genes and the so-called Clusters of Orthologous Groups (COGs). For the purpose of the analysis I am using the BioSAK COG2020 CLI but it looks like I have to align my DNA sequences differently and not by genes. Does somebody know how the sequences should be aligned for the CLI to produce meaningful results? The results set of COGs is empty. The problem is that these are nucleotide gene annotations from the Silva database. I suppose BioSAK expects annotated protein sequences but I have Nucleotide annotated sequences. The main question is: Is it possible to translate the Nucleotide to Protein sequences and annotate them to be used with BioSAK and do you know of a pipeline that does that?
Example command:
BioSAK COG2020 -m N -t 6 -db_dir ./COG_db_dir -i input.ffn

Example input Nucleotide fasta format:
>AB679109.1
GTGCCAGCAGCCGCGGTAATACGTAGGTGGCAAGCGTTGTCCGGAATTATTGGGCGTAA
AGCGCGCGCAGGCGGATCAGTCAGTCTGTCTTAAAAGTTCGGGGCTTAACCCCGTGATG
GGATGGAAACTGCTGATCTAGAGTATCGGAGAGGAAAGTGGAATTCCTAGTGTAGCGGT
GAAATGCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGACTTTCTGGACGAAAACT
GACGCTGAGGCGCGAAAGCCAGGGGAGCGAACGGGATTAGAAACCCCAGTAGTCC

Example from the COGs 20 fasta database:
>AAM04025_1 glutamate-1-semialdehyde 2,1-aminomutase [Methanosarcina acetivorans C2A]
MVSEVTLDKSRQMYEKAKTLIPGGVSSPVRAIKPYPFYTASADGSKIRDLDGNEYIDYCLAYGPAVLGHNHPVIKAAIKEQLDKGWLYGTPTELEVTLAEKVAGYYPSIDMLRFVSTGTEATMSALRLARGFTRKNKFIKIEGGFHGAHDAVLVKAGSGATTLGEPDSLGIPADFTKYTLQAPYNDIETMTTLVEKNRDDLAAVIIEPVLGNIGPILPLPGYLKELRKLTKENDVLLIFDEVITGFRLAMGGAQEYFGVVPDMTTLGKIVGGGLPIGVFGGRREIMEMIAPSGAVYQAGTFSGSPCSVAAGIAVLDYLKKEDIHAKLNSTGDYMRAVVSEIVEDEGLDYTVCGIASMFKIFFGAEPHNYQEALKCDKEGYLSFFHRMLANGVFLPPSQFETNFISAAHSEEDIEKTLEAYVENL

Please [edit] your question and give us a bit more detail. What is wrong with the output you are getting now? Do you get some sort of error message? What is it? Can you give us an example of you actual input (that's just a single fasta sequence, not an aligned file where we would need at least two). What exactly is your "aligned metagenome" file? What is it aligned against? Have you grouped by genes and aligned each gene? How did you identify the genes? Are they annotated or did you predict them? Are these bacteria?

Please show the first few lines of the results you're getting with your examples.

I have edited my question with more descriptive questions and info. The main question is can I extract COGs from Nucleotide-annotated sequences without additional protein sequencing?


text:
Upstream processing of raw single-cell sequencing data from bcl2fastq
I am trying to process some upstream data of single-cell sequencing, which are based on the 10x platform. They come directly from the output of bcl2fastq, so they are very raw and contain adapter sequences.
I first referred to the guide of cellranger (https://www.10xgenomics.com/support/software/cell-ranger/latest/analysis/running-pipelines/cr-gex-count). In the guide of cellranger, the steps to process these raw data seem to be simple. Basic steps are:
makefasta(bcl2fastq) --> cellranger count --> counts matrix --> other analysis

The questions I am concerned about are:

Is it necessary to apply standard NGS processing procedures, such as adapter removal and trimming of low-quality reads, which are not mentioned in the cellranger document?

I noticed that the R1 of fastq generated by the usual 10x genomic platform is mostly 26bp reads, which contains the information of umi index. But our data seem more complex, its looks have a ploy(T) tail and other sequence:

@A00199:1157:HKXC:2:1101:9914:1501 1:N:0:CACAATCCC+TTGTGGATA
ATCGTCCTCGGAGTGACGGTAATTCAGCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGGACTGGATTTTTATTTTTATTAAACTGGGGGTTTATTAAAAAAAGAAACAAGATGAACC
+
FFFFFFFFFFFFFFF:FFFF,:FFFFF::,FFF:F,F:FFF:F:FF,F:FFFFFFF,FF,,F,,,,,,,:FF,,,FFF,,,,,:,,,,,:,,,,,:F,:F,,,::,,,,,,,,,,,,,,,

Do I need to further process the R1 data?

Are there any classic related processes for reference?

Thanks for any suggestions.

Ask them https://community.illumina.com/s/

You are overthinking this.  Just use cellranger.  The stuff in R1 after the cell barcode and UMI is garbage. Cellranger knows that.

Thank you for your answer, so do I still need to do quality control like other NGS processes, such as adapter trim?

You don't really have to do adapter trim at all with RNASeq.  First of all, why would your reads be so short that you would read through to the adapter?  Second, STAR is a gapped aligner, it has no problem soft-clipping reads to get them to align.  Third, even if it didn't, what are the odds that it would matter?  You are just gene counting, so as long as a read gets counted as belonging to the right gene, it doesn't matter if there's cruft at the end.

Thank you for your answer, I belive you have answered my question very well. In fact, I just got involved in the analysis of scRNA-seq (I was previously engaged in genome analysis). We have completed the scRNA-seq from the sequencing company. The guide they provided suggesting checking and removing the adapter before downstream analysis, which made me confused because in my impression the insert size of the 10x platform is long than typic NGS. But I guess you have answered my question very well, thank you again.

If adapter sequences don't appear in the reference database, or won't be mapped to any sequences in the reference database, then removal is unnecessary for any quantification.
Adapter removal is also rarely needed prior to downstream processing with analysis suites, because they usually have any necessary filtering/removal built into the workflows.


text:
how to combine cpg information on both strands in analysis?
I’m trying to understand how bioinformatics tools import BS files, especially in terms of merging methylation data from the forward and reverse strands. How is this information represented in the final output? For instance, what happens when the CpG sites are symmetric, when data is available from only one strand, or when the two strands have differing methylation statuses? Thank you:).

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

What are BS files? Do you have any example files (e.g. the first few lines of a file)?

Hi, i have added an example pic :).


text:
Connecting and merging multiple RNA 3D structures in PyMol/Chimera
I have an optimized design for the corn RNA aptamer, PDB: 5BJO. I generated the 3D structures for the parts that I plan to add on trRosettaRNA. I want to connect them with the 5BJO PDB fiIe in PyMol, then perform Molecular dynamics. Is this possible?
This is the final design that I am trying to get. I have the 3D structure for each part separately. I was looking for a way to connect them together to build this model. Is it possible to assemble them in PyMol or other softwares? Thank you.

what do you mean by "Is it possible to connect and merge 4 different PDB files together," please clarify and produce an example of actual input and desidered output

@alpha do consider signing pippo1980's answer as accepted. Its good to do this all round.

@pippo1980
Thank you, I updated the description.

I see. Think you'll have to add manually the CONECT record between protein chains and RNA polymer. See Creating "stable" bond with PYMOL and the PDB file for 1ZR2:
Structure of a Synaptic gamma-delta Resolvase Tetramer Covalently Linked to two Cleaved DNAs

about pymol, when saving ALL (or ENABLED) objects in a .pdb file toghether it does keep atom numbering and remember that other Ligands (HETATM) are saved next to the chain they belong to

Pdb use ter records  this is something new to me too, try to use it before saving if you wish to create single chains out of multiple

@alpha do consider signing pippo1980's answer as accepted. Its good to do this all round and its a comprehensive answer including the points in the discussion.

in PyMol use, from menu :
File ---> Export Molecule... ---> and before hitting Save
select the molecules you want to save toghether (in the same .pdb file)
with the Selection drop down list either all (all the loaded mols are saved) or enabled (you can enable the ones you want to be saved toghether)
Examples of Save menu :

For command line see Thread: [PyMOL] independent positioning of two proteins, in preparation for RosettaDock?

Thank you. This can save the 4 structures in one PDB file, but I wanted to first connect them together, to assemble the final structure, then save into one PDB file. Alignment is not the way to do it, as far as I understand. Is it possible in PyMol or other softwares?

what do you mean by "Is it possible to connect and merge 4 different PDB files together," please clarify and produce an example of actual input and desidered output –
pippo1980