text:
Seurat Dimplot with different clustering IDs
In Seurat metadata I have assigned cells to some cell types with different resolutions
I have added the cluster identities to the object via Idents:
 Idents(control) <- a$predicted.celltype.l1

I then added another clustering in the same way:
 Idents(control) <- a$predicted.celltype.l2

Dimplot will now show the new identities:

How I can ask Dimplot to plot the clustering I want for example if I want to show predicted.celltype.l1 again (something like changing the default)?

The thing to tweak is the group.by argument of the Dimplot function.
It looks like you've already tried this in your attempts at solving a different problem.
The main thing you need to be aware of is that each new clustering needs to be saved as a separate metadata variable in order to do this. As I mentioned in one of my previous answers, I like to use the [[X]] notation for this:
control[["predicted.celltype.l1"]] <- a$predicted.celltype.l1
control[["predicted.celltype.l2"]] <- a$predicted.celltype.l2

Then you should be able to use both in the group.by argument:
DimPlot(control, reduction = "umap", group.by = "predicted.celltype.l1")
DimPlot(control, reduction = "umap", group.by = "predicted.celltype.l2")


text:
What does this Seurat argument mean
I have extensively read about percent mito in Seurat but I got more and more confused
Let's say we want to keep cells with maximum 20% of their expressed genes are mitochondrial genes, so shall I set percent.mito < 20 or percent.mito < 0.2 or percent.mito < 80 ?

Checkout the first Seurat tutorial. To filter cells with >20% mitochondrial counts:
obj[["percent.mito"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
obj <- subset(obj, subset = percent.mito < 20)

I tried reading the tutorial and things are not really clear there. @Exhausted, when the tutorial says "we filter cells", they mean filtering them out i.e. removing them. So in the subset call, one needs to state the criteria that the data needs to pass. It's quite clear that the percent is 1-100, not 0-1 though.

source: Stack Exchange

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

sort: votes


text:
Filtering cells from values in metadata
I metadata slot of Seurat object I have mapping score of each cell to a reference PBMC data like
pbmc <- subset(object, subset = object@meta.data$mapping.score > 0.5)
Error in FetchData(object = object, vars = unique(x = expr.char[vars.use]),  :
  None of the requested variables were found:
>
>
>

I want to only keep cells with mapping score > 0.5
Any help please?

your_new_object <- subset(your_object, subset = mapping.score > 0.5)

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


text:
Merging two dataframes in R
I have a file with cluster number of a scRNA-seq and corresponding annotated cell type like
> dput(head(anno))
structure(list(cluster = 1:6, celltype = c("T:CD4+NAIVE", "T:CD4+NAIVE",
"T:CD8+NAIVE", "NK:CD56+16+3+NKT", "T:CD8+NAIVE", "B:")), row.names = c(NA,
6L), class = "data.frame")

> unique(anno$cluster)
 [1]  1  2  3  4  5  6  8  9 11 12 15 16 19 20 21 22 23 24 25 27 28 29 30 31 35 36 37 38 39 40 41 42 43 44 45 46
[37] 47 48 49 50 51 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 80 81 82 83 84
[73] 85 86 87 88 89 90 91 92 93 94 95
> unique(anno$celltype)
 [1] "T:CD4+NAIVE"        "T:CD8+NAIVE"        "NK:CD56+16+3+NKT"   "B:"                 "NKT"
 [6] "T:CD8+EM"           "T:CD4+CM"           "mDC:"               "MACROPHAGE:"        "MONOCYTE:precursor"
[11] "NEUTROPHIL:"        "NK:CD56+16+3-"      "MAST:"              "MONOCYTE:"          "pDC:"
[16] "T:Reg"

and another data frame with cell IDs and corresponding cluster number of each cell
> dput(head(df2))
structure(list(cluster = c(1L, 1L, 1L, 1L, 1L, 1L), cellid = c("cell17203",
"cell17205", "cell17206", "cell17207", "cell17208", "cell17209"
)), row.names = 15498:15503, class = "data.frame")

> unique(df2$cluster)
 [1]  1  2  3  4  5  6  8 12 16 19 21 23 27 29 30 35 36 40 53 55 92 11 24 38 41 46 90  9 15 20 47 74 57 39 42 25
[37] 56 64
>
> dim(df2)
[1] 16056     2
>

I want to create a cell type column in the second data frame (df2) in which relate cell type from the first data frame to cluster number for instance cluster 1 is T:CD4+NAIVE
I have tried these
> df2[anno, on = "cluster"]
Error in `[.data.frame`(df2, anno, on = "cluster") :
  unused argument (on = "cluster")
> df2[c, on = anno$cluster]
Error in `[.data.frame`(df2, anno, on = c$cluster) :
  unused argument (on = anno$cluster)
> df2[c, on = df2$cluster]
Error in `[.data.frame`(df2, anno, on = df2$cluster) :
  unused argument (on = df2$cluster)

You used data.table notation without actually having data.table:
library(data.table)
setDT(anno)
setDT(df2)
df2[anno, on = "cluster"]

    cluster    cellid         celltype
 1:       1 cell17203      T:CD4+NAIVE
 2:       1 cell17205      T:CD4+NAIVE
 3:       1 cell17206      T:CD4+NAIVE
 4:       1 cell17207      T:CD4+NAIVE
 5:       1 cell17208      T:CD4+NAIVE
 6:       1 cell17209      T:CD4+NAIVE
 7:       2      <NA>      T:CD4+NAIVE
 8:       3      <NA>      T:CD8+NAIVE
 9:       4      <NA> NK:CD56+16+3+NKT
10:       5      <NA>      T:CD8+NAIVE
11:       6      <NA>               B:

For merging dataframes, I find it easiest to use the tidyverse / dplyr functions inner/full/left/right_join. See the "Data Transformation Cheatsheet" on this page. For the merge that @user438383 has mentioned, this would be a left_join:
library(tidyverse)
anno <- data.frame(cluster = 1:6, celltype = c("T:CD4+NAIVE", "T:CD4+NAIVE",
                   "T:CD8+NAIVE", "NK:CD56+16+3+NKT", "T:CD8+NAIVE", "B:"))

df2 <- data.frame(cluster = c(1L, 1L, 1L, 1L, 1L, 1L),
                  cellid = c("cell17203", "cell17205", "cell17206",
                             "cell17207", "cell17208", "cell17209"),
                  row.names = 15498:15503)

> left_join(anno, df2)
Joining, by = "cluster"
   cluster         celltype    cellid
1        1      T:CD4+NAIVE cell17203
2        1      T:CD4+NAIVE cell17205
3        1      T:CD4+NAIVE cell17206
4        1      T:CD4+NAIVE cell17207
5        1      T:CD4+NAIVE cell17208
6        1      T:CD4+NAIVE cell17209
7        2      T:CD4+NAIVE      <NA>
8        3      T:CD8+NAIVE      <NA>
9        4 NK:CD56+16+3+NKT      <NA>
10       5      T:CD8+NAIVE      <NA>
11       6               B:      <NA>


text:
Error when trying to move multiple file types in wdl
I tried to add a "gather_results" task to a workflow that takes in an array of Files and moves output files to a new folder of a specified name.
Cannot interpolate Array[File?]+ into a command string with attribute set wdl PlaceholderAttributeSet(None,None,None,Some( ))]
The wdl file passes womtools and I have no idea what a PlaceholderAttributeSet is, which makes debugging quite difficult.
--
gather files task
task gather_result{
  input{
    String dir
    String? subDir
    Array[File?]+ files
  }

  command
  {
    mkdir --parents ${dir}/${default="" subDir}; /bin/mv -t ${dir}/${default="" subDir} ${sep=" " files}
  }
  runtime
  {
    docker_image: "pipeline"
    memory: "2 GB"
  }
  output
  {
    String out = stdout()
  }
}

How it's called:
  call gather_result as gather_files {
      input: files=[create_hic.output_hic[1],
                    dedup.deduped_bam,
                    calculate_stats.stats[0],
                    calculate_stats.stats_hists[0],
                    calculate_stats.stats_json[0]],
            dir=OutputDir
    }

WDL can be a pain. Try and simplify it (reduce input complexity by adding one piece of files at a time). Wil need trial and error. Drop the ? if U can. Check all input types. One might need a flatten.

If I drop the ? I get a womtools error saying wld can't interpolate Array[File]+ to Array[File?]+

I think you'll need to convert your Array[File?]+ to a Array[File] before you can use ~{sep=" " files} in the command section.
You can use the select_all function for that, which will only keep valid files.
https://github.com/openwdl/wdl/blob/main/versions/development/SPEC.md#arrayx-select_allarrayx


text:
Low pass sequencing has been reported to detect common variants. How low can one go and get reliable data? Is 2X pass sequencing analysis possible?
I would like to use low pass sequencing to replace a genotyping chip to be able to detect variants up to 0.1 % allele frequency in available population data. What is the minimum depth I can opt for to retain the confidence of base calls? Can GATK be used to analyse this data?

This is a good idea - low pass sequencing has shown to be better than genotype arrays for things such as GWAS or QTL mapping.
The primary way of processing low-pass sequencing data is to use imputation - partciularly methods such as GLIMPSE or QUILT and a large reference panel such as the HRC or TopMed in order to refine the genotype calls.
This is a figure from the GLIMPSE paper showing how you can retain variants of different frequencies at different levels of sequencing depth.

And a similar image from the QUILT paper:

You might also want to check out some of the methodology in this paper http://m.genome.cshlp.org/content/31/4/529


text:
How to quantifiy of specific genes from shotgun metagenome?
I have googled a "lot", couldn't find any specific answer to the question. So, I am here seeking for your guidance. My question is similar to this. I have several metagenome (n=30). But for simplicity’s sake, let's assume only one metagenome. I wanted to quantify lasA genes in the metagenome. So far, what I have done (not sure if it's the right approach) is:

Download the lasA genes from uniport.
I then made the multiple sequence alignment (lasA.msa) of the lasA gene downloaded in the previous step.
The protein profile was generated hmmbuild res.hmm lasA.msa
Genome assembly (contigs) were done with megahit megahit -1 r1.fastq.gz -2 r2.fastq.gz --min-contig-len 1000 -m 0.85 -o ASSEMBLY/ -t 20
ORFs were predicted and extracted using getorf from the contigs (from step 4)
The profile (from step 3) was used to search the ORFs using hmmsearch
Several hits were obtained from hmmsearch indicating that the lasA gene is present in the metagenome

Now, my question is, how do I count the number of reads affiliated to the lasA gene? I feel that I need to map the reads to the assembly, but I don't know how to go about that.
NOTE: I have used MG-RAST, Humann3, and eggNOG-mapper, but none of them showed the presence of lasA gene, however the above-mentioned hmmsearch was able to find the gene, what could be the reason? Problem with my approach? Or hmmsearch using the protein profile is more sensitive?
Thanks and regards

The way I would do it is to make an assumption that the reads mapping to your contig == (abundance of the genes in the contig)/contig length.

Map raw reads to the contig containing your hit (original contig, pre-translation)

Normalise it with RPKM:

RPKM = NumberOfReadsThatAligned / ( ContigLength/1000 * totalNumberReads/1,000,000 )

This will give you the relative abundance of the contig which is an APPROXIMATION of the abundance of your gene.

Alternatively, to be more precise you can FragGeneScan your raw reads, extract the sequence of your putative protein and map translated raw reads to it.
For your other question: Yes, HMMs are more sensitive than blast or diamond since its output is probability score that these protein contain motifs extracted from your alignments (given that there is one). For that reason you need to look at the bitscore and E-value. If the bitscore is low and/or E-value is high, you have a high chance of false-positive.
Best wishes,
A


text:
Python script to calculate properties of amino acids and storing them in matrix
I want to make a python script that can get as input sequence of N amino acids, and the output is a matrix of Nx6 containing 6 features.
The features are:

Computed volume
Hydrophobicity
Surface charge
Secondary structure (SS)
Go, InterPro and Pfam type

I have developed this script :
#Importing Libraries
import numpy as np
import pandas as pd
from Bio.SeqUtils.ProtParam import ProteinAnalysis
from Bio.SeqUtils.ProtParam import ProtParamData
from quantiprot.metrics.aaindex import get_aa2volume, get_aa2hydropathy
from quantiprot.metrics.basic import average

#Input amino acid
 my_seq = input("Enter Your Sequence : ")
 SEQ = my_seq.upper()

#Analysis
analysed_seq = ProteinAnalysis(SEQ)
MW=  analysed_seq.molecular_weight()
Gravity = analysed_seq.gravy()
aa_composition = analysed_seq.count_amino_acids()
aa_percentage=analysed_seq.get_amino_acids_percent()
HP = analysed_seq.protein_scale(window=7, param_dict=ProtParamData.kd)
Sec_Str = analysed_seq.secondary_structure_fraction()

 #Amino Acid Compsoition Calculation
  print("\n" ,"AA count:",aa_composition,"\n")

  #Molecular Weight Calculation
   print("Molecular Weight : ",MW)

   #Gravity Calculation
    print("\n","Gravity:","\n\n",Gravity)
    #Hydrophobicty Calculation by using kd scale
    print("\n","Hydrophobicity:","\n\n",HP)
    #kd → Kyte & Doolittle Index of Hydrophobicity
   #Flex → Normalized average flexibility parameters (B-values)
   #hw → Hopp & Wood Index of Hydrophilicity
   #em → Emini Surface fractional probability

  #Amino Acid Percent Calculation

 print("\n","Percentage of Amino Acids in Protein:","\n\n",aa_composition)

 #Volume Calculation by using quantiprot library
 vol =  get_aa2volume(analysed_seq)
print("\n","Volume of amino acids in protein:","\n\n",str(vol))

I have calculated volume, Hydrophobicity, and SS but I don't know how can I find the

Go, InterPro and Pfam type
polarity, i.e. surface charge

Type and polarity are static, right? So just import from some db. Surface accessibility changes based on folding doesn't it? So unless U do something really fancy with alphafold that might be hard...

RSA can be computed with NACCESS which I think also has a BioPython interface

RSA requires a structure. Predicted SS does not, but assigned SS does and is  is folds better (bad pun intended). You have also Normalized average flexibility parameters (B-values), which requires both a structure and metadata on that structure as the B-factors depend on the resolution (unless you have a protocol for it, it is not a rabbit hole worth going down, not due to the maths, but due the isotropic and anisotropic business).

Which dB I can use to import type and polarity?

If by "Polarity", you mean global charge, then that's isoelectric point (pI), which is actually useful for the protein engeering folk (pI close to pH is less soluble). If you mean some sort of per-residue value based on charged and polar residues in a window, then that would be something nearly similar to a Kyte & Doolittle Index. The per-voxel surface charge on a protein, requires a structure (and the Advance Poison-Boltzmann Solver).

how can I find the type?

Type is vague word. Take GNB1, the G-protein subunit beta-1. It has a GO term GO:0003924 (GTPase activity), an InterPro  IPR001632 (G-Protein B), and PFam 7x PF00400 (WD40 repeats). So the word describing what type of protein a protein is depends on the classifier in question.

but I want to find the type if any input sequence?

Read Matteo's comment - "type" is a vague word. Define it better; give it more context and maybe someone will be able to help you.


text:
How are tertiary protein structures represented in computers?
What data structure is used in representing protein structures in computers so that we can apply algorithms?
Matrix or Graph or tree?
P.S. I am absolutely new to Bioinformatics.

the same was asked on SO: https://stackoverflow.com/questions/68452847/how-are-tertiary-protein-structures-represented-in-computers

This is too vague to answer.  A small compound can be represented with a string (SMILES, no 3D data though), a matrix (mol file) or a graph network object (eg. a RDKit Chem object), yet some methods return matrices (e.g. distance). The same goes for a protein...

@MatteoFerla, what more info do u need to make this question clearer?

You kind of have answered the question in the questions. Both matrix and graph are used to represent structures in memory, but they have different uses. Some algorithms will use one or the other.
A PDB, CIF, etc. file is basically a matrix presentation.
But a major pressure for a graph network is if they do not operate in cartesian space (x,y,z) but torsion/dihedral space (d,φ,ψ).

comments from above fleshed out.

Both: A graph network is a handy way to represent a molecule, but a matrix makes for fast maths.
A graph network is a list of nodes (atoms), a list of edges, the connections between nodes and a list of properties associated with the nodes (symbol, name, index, residue index, residue number, chain, coordinates etc.) and occasionally a list of those of the edges (single, double etc.).
A protein is a big molecule. Same rules apply.
A small compound can be represented with a string, i.e. SMILES, which is a form of graph network even if it looks nothing like one.

In a mol file and its sdf/mol2 derivatives, a tabular representation is used to store coordinates (properties) of the atoms of the small molecule (top part below), while the connectivity has an edge-node relationship, i.e. are edges. Note that the index of the properties in the first block is the id of the atom.
     RDKit          3D

 12 12  0  0  0  0  0  0  0  0999 V2000
   -5.5310   24.4560   -8.2060 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.9560   21.4190   -7.6300 O   0  0  0  0  0  0  0  0  0  0  0  0
   -5.4170   23.9250   -6.8750 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.2770   19.3290   -8.7330 O   0  0  0  0  0  0  0  0  0  0  0  0
   -3.9460   23.5670   -6.5510 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.3270   19.7590   -6.4400 O   0  0  0  0  0  0  0  0  0  0  0  0
   -3.4590   22.5590   -7.5870 N   0  0  0  0  0  0  0  0  0  0  0  0
   -3.5890   23.1360   -8.9920 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.0600   23.5190   -9.2240 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.0390   22.1150   -7.2780 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6090   21.0030   -8.2230 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0430   20.3530   -7.7350 S   0  0  1  0  0  0  0  0  0  0  0  0
  1  3  1  0
  1  9  1  0
  3  5  1  0
  9  8  1  0
 12  2  1  1
 12  4  1  0
 12  6  1  0
 12 11  1  0
  5  7  1  0
  7  8  1  0
  7 10  1  0
 10 11  1  0
M  END

A PDB, CIF, etc. file is the same. It has a per-atom tabular representation of the properties, with an optional connectivity part (CONECT entries) for non-standard residues (not AAs and DNA).
Now when these are loaded into memory, they become an array of atoms and bonds, each with their properties —graph network.
However, the coordinates within get often converted into a matrix for speedy operations.
For example to calculate an atom to atom distance matrix, converting the coordinates with the array of atoms to a matrix can be very easily be converted into a distance matrix via simple matrix calculations.
So both matrix representations and graph network representations are used.
For 3D visualisation, a pre-existing library is used and bonds and atoms become primitive shapes, whose splines may be connected to each other (in knots/vertices) or simply the bonds shapes clip through the atoms. That is there is not graph relationship.


text:
How to download the whole genome of an organism from NCBI using R or Python?
I need to download some of the genomes of human gut microbiome in FASTA format. Anybody has any script to download?

I'd suggest the bit-dl-ncbi-assemblies from the Bioinformatics Tools package.
$ bit-dl-ncbi-assemblies

 --------------------------------  HELP INFO  ---------------------------------

  This program downloads assembly files for NCBI genomes. It takes as input
  assembly accessions (either GCA_* or GCF_*) and optionally a specification of
  which format to download (currently: genbank, fasta, protein, or gff). For
  version info, run `bit-version`.

    Required input:

        - [-w <file>] single-column file of NCBI assembly accessions

    Optional arguments include:

        - [-f <str>] default: genbank
                  Specify the desired format. Available options currently
                  include: genbank, fasta, protein, gff, feature_tab, report, stats.

        - [-j <int> ] default: 1
                  The number of downloads you'd like to run in parallel. Write speeds
                  can become problematic generally with around 15 or more, so a max
                  of 12 will be used to be safe even if more are requested.

    Example usage:

     bit-dl-ncbi-assemblies -w ncbi_accessions.txt -f protein -j 4

Alternatively, you can use this code from the Biopython cookbook:
import os
from Bio import SeqIO
from Bio import Entrez

Entrez.email = "A.N.Other@example.com"  # Always tell NCBI who you are

seq_ids = ["EU490707"] # Your accessions go in this list
for sid in seq_ids:
    filename = f"{sid}.fasta"
    if not os.path.isfile(filename):
        # Downloading...
        net_handle = Entrez.efetch(
            db="nucleotide", id=sid, rettype="fasta", retmode="text"
        )
        out_handle = open(filename, "w")
        out_handle.write(net_handle.read())
        out_handle.close()
        net_handle.close()
        print("Saved")


text:
Doing functional enrichment of microbiome data
I have done gene set enrichment analysis which is now straight forward. Such as taking gene set and making a query. For microbiome data it is not that straight forward which I found out.
Step done so far..

Ran Kraken2
Generated OTU table
Downstream analysis using phyloseq.

Now if I have to the functional enrichment from start I would need to use qiime2 which I didn't do earlier.
Now the other option after reading through paper I found kegg can also perform the microbiome functional enrichment
ISSUE
Unlike other database I'm not sure where do I put query in KEGG for the bacterial genus level whatever found.
Any suggestion how to proceed would be really helpful


text:
DESeq2. Which is better: using the altHypothesis argument in the results function or manually filtering for your desired P-value and Log2FoldChange?
In DESeq2 you can identify your differentially expressed (DE) genes using the results function. I noticed people in my lab using the results() function with the minimum number of arguments supplied (below) and then converting the resulting object to a data frame filtering for the genes they considered DE (where they considered a fold change of +/- 2 and a p-value of 0.05 to be DE).
res1 <- results(dds, contrast(condition, treated, untreated))

Alternatively, I noticed that the results function allows you to set both the P-value and minimum fold change when you use the 'altHypothesis' argument (below). The resulting list from using this method is considerably shorter than the list from manually filtering for the same P-value and fold change. I am aware that that by using the 'altHypothesis' argument it comparisons are made using a Wald test rather than a LRT test, however, I am completely unfamiliar with these two tests and I was not able to understand their significant differences in the context of RNA-seq analysis.
res1 <- results(dds, contrast(condition, treated, untreated), alpha = 0.1, lfcThreshold = 1, altHypothesis = "greaterAbs")

Due to the difference in size of the lists, and the goal of being consistent moving forward for the lab, does anyone have any insight as too which one is better? Neither myself or my lab are particularly used to this analysis but we would like to try to follow best practices as much as possible.
Thank you

I love that Deseq2 has altHypothesis="lessAbs" !!!!!!! I've used it a ton for RNA-seq. However, now I'm working with mass spect data using limma. Is there a way to make limma do the same thing?

@MaryAllen welcome to the site. You might consider asking the question in your comment as a separate question. There is a fair bit of popularity for DSeq2 here

If you have an existing understanding of a fold change that you would consider to be biologically relevant, it will be better to use the lfcThreshold for comparisons. Working out that fold change is the hard bit.
When using the Wald test, the p-value is associated specifically with that threshold, rather than the difference from zero.
This is explained in the DESeq2 Analysis Vignette:

It is also possible to provide thresholds for constructing Wald tests of significance. Two arguments to the results function allow for threshold-based Wald tests: lfcThreshold, which takes a numeric of a non-negative threshold value, and altHypothesis, which specifies the kind of test. Note that the alternative hypothesis is specified by the user, i.e. those genes which the user is interested in finding, and the test provides $p$ values for the null hypothesis, the complement of the set defined by the alternative. The altHypothesis argument can take one of the following four values, where $β$ is the log2 fold change specified by the name argument, and $x$ is the lfcThreshold.

greaterAbs - $|β|>x$ - tests are two-tailed
lessAbs - $|β|<x$ - p values are the maximum of the upper and lower tests
greater - $β>x$
less - $β<−x$


text:
How to download data from SRA in Linux systems via the command line?
My workflow for downloading data from SRA has been the following:

Access SRA Run Selector.
Enter the accession number for the project of interest.
Download "Accession List" for the "Total" data. The accession list is basically, as the name says, a .txt file listing all accession numbers for each run in a given bioproject.
For each accession number in the accession list I run the following command:

$ fastq-dump --accession <accesion_number> --outdir .

How to download data from SRA in Linux systems via the command line? What are fast and scalable ways of achieving this? What should one know before downloading data?
Although there are some similar discussions in this StackExchange (eg. Trying to retrieve sra data using sra explorer, Downloading dataset from SRA (SOLiD Platform), Downloading SRA Files from AWS, etc).
The discussion on What's the best way to download data from the SRA? Is it really this slow?
 seems like a good starting point, but I wasn't able to find a straightforward guide on how to download data from SRA, e believe that a well thought and complete answer would help beginners alike. I am preparing an answer to the problem myself, but I thought that the community would probably give some interesting insights (:

I wouldn't say the download itself is slow but in my experience sra to fastq conversion took much more time than the download, and using fasterq-dump instead of fastq-dump did not help much. I would take a look at ENA (https://www.ebi.ac.uk/ena/browser/home) and see if they offer fastq files for the dataset(s) you are interested in.

The dataset I am looking for is actually available at ENA, and indeed, the workflow is way nicer using their interface.

So far, grabseqs is the easiest option, a wrapper for fastq_dump and fasterq-dump, you can install it with conda, I recommend you to use an environment:
conda create -n grabseqs -c louiejtaylor -c bioconda -c conda-forge grabseqs

conda activate grabseqs

grabseqs sra SRRNNNNNNNN

conda deactivate

And download directly projects, runs, and bioprojects from SRA with grabseqs sra SRRNNNNNNNN. Check it out here: https://github.com/louiejtaylor/grabseqs
You can install in a different way just in case the first one doesn't work:
conda create -n grabseqs -yc bioconda parallel-fastq-dump pigz
conda activate grabseqs
pip install grabseqs

grabseqs sra SRRNNNNNNNN

conda deactivate

Nice, I will test this in a couple of hours ^^

I didn't manage to install it easily, it is presenting some kind of conflict between the dependencies. I will try to find something else.

I did a slight modification just in case the official installation steps didn't work.


text:
Converting Fastq files to Fasta files on Ubuntu
I am new to bioinformatics and programming. I tried to convert my fastq files to fasta file, but got this. What could be wrong?
(base) jszheng@Greater:~/seqtk$ seqtk seq -a KF0124_1_1.fastq > KF0124_1_1.fasta
[E::stk_seq] failed to open the input file/stream.

Answer from @haci, converted from comment:
You should run your command within that directory (the directory where the fastq file is), this does not seem to be the case as far as I see. Alternatively you can provide the absolute path of your fastq file.
cd into the directory where the fastq file is (you should be able to see your file of interest when you type ls) and then run the command in your question:
cd ./genome_diirt
seqtk seq -a KF0124_1_1.fastq > KF0124_1_1.fasta

Alternatively you can provide the absolute path of your fastq file:
seqtk seq -a ./genome_diirt/KF0124_1_1.fastq > ./genome_diirt/KF0124_1_1.fasta


text:
programmatically select and label peptides within pdb file
I have a pdb file of a protein structure, and I would like to programmatically select certain peptides and label them in an annotated pdb file output, so that when the output pdb file is opened again by one of the 3D viewers, one can highlight the annotated peptides.
I can do this interactively using UCSF Chimera or https://www.rcsb.org/3d-view,
and with Chimera, I can even save a session which is a python file that contains the lines of my manual annotations as lists of coordinates.
whereas in RCSB, saving a session produces a molx compressed file that includes hashes such as this one, where myseq is defined by a range within the file:
              {
                "parent": "yNaD8QxiUNHaHCvxrC5rew",
                "transformer": "ms-plugin.structure-component",
                "params": {
                  "type": {
                    "name": "bundle",
                    "params": {
                      "hash": 1673297064,
                      "elements": [
                        {
                          "groupedUnits": [
                            [
                              0
                            ]
                          ],
                          "set": [],
                          "ranges": [
                            1669,
                            1777
                          ]
                        }
                      ]
                    }
                  },
                  "nullIfEmpty": true,
                  "label": "myseq2"
                },

So I want to do the same programmatically, e.g.:
annotate_sequences --input file.pdb --sequences 1:ARNDC,2:CDNRA,foo:VYWTSP

I expect an output file where these three peptides are annotated with the labels 1, 2, and foo, so when someone opens the file in a 3D structure viewer, then can select/deselect the annotations.
Any ideas on how to do this?

Do you mean something specific by annotated pdb file? Regular pdb files don't have place for annotations of chains (peptides).

Not as REMARK? Is there another format for this?

There is no consensus on how different programs do representations (example of NGL vs. PyMOL) or what they are called so there is not official way to store selections.
Annotations can be REMARK lines, which have a specific format or simply lines after a hash at the bottom of the file get ignored by most viewers but that specific program may read. Here is a scored Rosetta PDB file for example:
HEADER                                            14-JUN-21   XXXX
EXPDTA    THEORETICAL MODEL
REMARK 220
REMARK 220 EXPERIMENTAL DETAILS
REMARK 220  EXPERIMENT TYPE                : THEORETICAL MODELLING
REMARK 220  DATE OF DATA COLLECTION        : 14-JUN-21
REMARK 220
REMARK 220 REMARK: MODEL GENERATED BY ROSETTA
REMARK 220  VERSION 2020.49+release.fac69dd
ATOM      1  N   MET A   1     107.793 134.398   6.022  1.00  0.00           N
ATOM      2  CA  MET A   1     108.468 134.162   7.299  1.00  0.00           C
ATOM      3  C   MET A   1     108.644 135.405   8.219  1.00  0.00           C
ATOM      4  O   MET A   1     108.050 135.416   9.303  1.00  0.00           O
ATOM      5  CB  MET A   1     109.782 133.398   7.080  1.00  0.00           C
ATOM      6  CG  MET A   1     110.650 133.262   8.335  1.00  0.00           C
ATOM      7  SD  MET A   1     112.224 132.454   8.016  1.00  0.00           S
...
# All scores below are weighted scores, not raw scores.
#BEGIN_POSE_ENERGIES_TABLE
label fa_atr fa_rep fa_sol fa_intra_rep fa_intra_sol_xover4 lk_ball_wtd fa_elec pro_close hbond_sr_bb hbond_lr_bb hbond_bb_sc hbond_sc dslf_fa13 omega fa_dun p_aa_pp yhh_planarity ref rama_prepro total
weights 1 0.55 1 0.005 1 1 1 1.25 1 1 1 1 1.25 0.4 0.7 0.6 0.625 1 0.45 NA
pose -17624.7 2135.24 10828.2 38.8161 591.762 -521.596 -5634.4 9.63107 -1337.76 -585.258 -616.219 -453.443 0 197.59 3291.97 -602.3 2.74143 1088.19 99.1601 -9092.4
MET:NtermProteinFull_1 -2.63316 0.99093 2.70981 0.01206 0.04738 -0.17936 0.08728 0.00151 0 0 0 0 0 0.35192 1.9909 0 0 1.65735 0 5.0366

Hacked b-factors
Another way is to store numerical data as b-factors. The AlphaFold2 data for example stores pLDDT in the b-factor column, while consurf stores conservation score there. So a boolean 1=show, 0=hide could be done for example, which would allow one to select these, e.g. in PyMOL:
# store chosen selection...
alter chosen, "b=1"
sort
save foo.pdb

# load it back...
select chosen, b<0

Although this makes no sense as one can save a session file...
Also note, that for multiple selections, you'd need create different models, but whatever happens the other non-coordinate column, occupancy, should not be touched.
Sharing
If the problem is that you want to share a protein with someone and want to let the recipient know about features of interest but are using different software than you, you may want to consider Michelanglo, a website I wrote that allows one to share a webpage with an interactive protein viewer (using NGL.js, hence the name) and a side description with custom annotations that control the protein representation (example).

How would one annotate the certain peptides within the sequence the PDB corresponds to? E.g. starting with the full sequence + PDB structure, how do we end up with an annotated structure for peptides within, that is then shared via Michelanglo?

The hacky B-factor annotations I mentioned have nothing to do with Michelanglo, which has a GUI to add protein-controlling links in the description of a page a registered user has created by uploading a regular structure or gene name —i.e. after upload.


text:
Where can I find a basic example for testing association between a phenotype and each of the variants of a chromosome?
I am a statistics student trying to find data suitable for a particular method. From what I have read so far, it is of interest to determine if variants within a chromosome are associated with a phenotype. I would like to understand this more and if possible, find an example of this application. In particular, what test statistic is used when associating each chromosome variant to a phenotype? Also, is it possible to gain access to examples of SNP Chip data fitting this context?

There are a very large number of such methods. I would recommend googling "GWAS" or "genome wide association analysis" and you will find many many results.
One relatively simple example would be this plant genotype-phenotype dataset, associated with this publication. Some basic summary description of the GWAS dataset can be found here. By now the methods in this paper are very old and naive seeming, but it gives an idea for how people approached this problem in the golden age of GWAS. The data are small and simple and full of strong associations. Additionally, the phenotype measurements are based on isogenic inbred lines, which are more or less homozygous, making the analysis much easier, interpretation a lot simpler.
Another point in favor of using the plant data is that there is actually some strong genotype-phenotype association with a really strong mechanistic basis, like the various flowering time associations. For better or for worse, such unambiguous associations are often absent in e.g. human phenotypes of interest.
For some background on the general problem of controlling for population structure in GWAS, see the old EMMA paper. A much more recent example of current statistical methodology that has been making the rounds on Twitter can be found here.
For very large catalogs of GWAS data in humans, I can point you here or here. The UK Biobank dataset in particular is very commonly used and quite high quality.


text:
SCTransform Warning: in theta.ml(y = y, mu = fit$fitted) : iteration limit reached
Running SCTransform on my seurat object produces the warning:
Warning in theta.ml(y = y, mu = fit$fitted) : iteration limit reached

What does this warning mean? Can I trust the output?
SCTransform runs through the data producing this warning until the progress bar reaches 100%. The resulting seurat object does have an SCT assay and I have successfully done downstream analysis with it. I am concerned that my SCTransform has not successfully done what it is supposed to do. Is there a way to cross-check?

The explanation can be found in link posted by ATpoint in his comment:

I get at least 10 posts on GitHub, GitMemory and Biostars just on the top page.

https://github.com/satijalab/seurat/issues/1378
https://github.com/satijalab/seurat/issues/2578
https://github.com/satijalab/seurat/issues/1378
https://www.biostars.org/p/450394/

The first author of the SCTransform paper, Christoph Hafemeister had posted this answer in a similar query on the satijalab GitHub page https://github.com/satijalab/seurat/issues/1378
"These warnings are related to the estimation of the theta parameter of the Negative Binomial regression model used to describe the molecule count distribution of each gene. When a gene has very few non-zero observations estimation of that parameter is not reliable and that's what the message indicates. This is not a problem, since the estimated theta parameters are not used directly for normalization but are regularized (sharing information across genes and thus lowering uncertainty)."


text:
Help with Seurat QC ambiguity
I have four PBMC samples from 10X scRNA-seq
> cancer
An object of class Seurat
36601 features across 18338 samples within 1 assay
Active assay: RNA (36601 features, 0 variable features)
>

And this is Seurat QC plot

If I am not wrong, samples were of relatively low quality, with gene expression data revealing the presence of mitochondrial genes, as well as MALAT-1, which are suggestive of poor sample quality (dead/dying cells).
Any way by this plot I filtered cells to remove data that had more than 20 "percent" mitochondrial expression , > 2000 features or < 100 features
cancer <- subset(cancer, subset = nFeature_RNA > 100 & nFeature_RNA < 2000 & percent.mt < 20)

By this I lost most of cells
> cancer
    An object of class Seurat
    36601 features across 6,883 samples within 1 assay
    Active assay: RNA (36601 features, 0 variable features)
    >

Please, can somebody experienced in scRNA-seq inspect this plot and my thresholds and tell me if I am wrong, any suggestion, if losing this number of cells is normal
Thank you

Answer from @haci, converted from comment:
I think you should not be using a cut-off higher than 20 or 25 for the reasons I explained above.
I would advise against using cells with more than 20-25% percent.mito. This metric is used as a proxy for cell death and there is not much point keeping these cells if the cell death is so severe to the point where let's say half the reads map to the mitochondrial genome (with a 50% cut-off). There will be many caveats one being apoptotic cells from different lineages clustering together, not interesting unless you are studying apoptosis related phenomena. In my opinion 6000 good quality cells are better than 10000 apoptotic cells.

From the graphs, you're losing most of the cells through the percent.mitochondrial threshold. The distribution is so broad on that mitochondrial plot that I probably wouldn't filter for mitochondrial expression (except possibly at ~95%). If I did have a lower threshold, I'd put it at around 50% rather than 20%, because that's where there's a change in density on the mitochondrial plot.
This suggests that there's a very large mitochondria-rich population in the PBMC dataset. I don't have a good understanding of biology and cell populations, but if I were presented with these data I'd probably try clipping out those cells with 65-95% mtDNA and see if there is a strong enrichment for markers that correspond to any particular cell population (e.g. platelets).


text:
Parallelizing microRNA targets
I  am trying to look for miRNA targets using a file called Conservedfamily.txt from the Zebrafish target scan fish website. I have written a python program to extract each of the miRNAs which works fine. However, for a big list of microRNAs, I must provide a name every time in line 22 and get one individual output file. I am thinking is there any option to parallelize the process so that somehow, I can provide a list of all microRNAs and obtain all my targets at the same time. The python code is attached below. I have posted the question in biostars but go no reply from them for 4 days. Hence I am posting the question in Stack exchange. The python code is attached below.
#!/usr/bin/python

#Import relevant modules

import sys

conservedfamily = open(sys.argv[1], 'r')
targetgenelist = open(sys.argv[2], 'w')
targetgeneid = open(sys.argv[3], 'w')

geneid = []

targetgenelist.write('miR family' + '\t' + 'Gene ID' + '\t' + 'Gene Symbol' + '\t' + 'Transcript ID' + '\t' + 'Species ID' + '\t' + 'UTR start' + '\t' + 'UTR end' + '\t' + 'MSA start' + '\t' + 'MSA end' + '\t' + 'Seed match' + '\n')

for line in conservedfamily:
  splitline = line.strip().split('\t')

  if splitline[0] == 'miR-210':
    targetgenelist.write(line)
    if not splitline[1] in geneid:
      geneid.append(splitline[1])
      targetgeneid.write(splitline[1] + '\n')

conservedfamily.close()
targetgenelist.close()
targetgeneid.close()

python miR99target.py 'Conserved_Family_Info.txt' 'miR-210targetlist.txt' 'miR-210targetid.txt'

I'd suggest loading in all of the miRNA tags as a set(). Then, you can just check the set membership each time.
#!/usr/bin/python

#Import relevant modules

import sys

conservedfamily = open(sys.argv[1], 'r')
targetgenelist = open(sys.argv[2], 'w')
targetgeneid = open(sys.argv[3], 'w')

# Create a set of all miRNA tags. Each tag should be on its own line.
with open(sys.argv[4], 'r') as miRNAlist:
    miRNA_set = {line.strip() for line in miRNAlist}

geneid = []

targetgenelist.write('miR family' + '\t' + 'Gene ID' + '\t' + 'Gene Symbol' + '\t' + 'Transcript ID' + '\t' + 'Species ID' + '\t' + 'UTR start' + '\t' + 'UTR end' + '\t' + 'MSA start' + '\t' + 'MSA end' + '\t' + 'Seed match' + '\n')

for line in conservedfamily:
  splitline = line.strip().split('\t')

  # Check if current tag is in the desired set
  if splitline[0] in miRNA_set:
    targetgenelist.write(line)
    if not splitline[1] in geneid:
      geneid.append(splitline[1])
      targetgeneid.write(splitline[1] + '\n')

conservedfamily.close()
targetgenelist.close()
targetgeneid.close()

One additional comment:
I'd suggest using the python module for csv reading/writing, even for quick and dirty scripts as well as using context managers for reading/writing files.


text:
Seurat heatmap for two conditions
I have a Seurat object of four cancers and four controls
> colnames(pbmc@meta.data)
 [1] "orig.ident"       "nCount_RNA"       "nFeature_RNA"     "samples"
 [5] "Condition"        "Gender"           "percent.mt"       "nCount_SCT"
 [9] "nFeature_SCT"     "SCT_snn_res.0.75" "seurat_clusters"  "SCT_snn_res.0.5"
> unique(pbmc@meta.data$Condition)
[1] "control" "cancer"
>

I want to have a heat map of differentially expressed between cancers and controls like figure 3A in this paper
https://www.frontiersin.org/articles/10.3389/fgene.2020.610682/full

Do you have any idea please how I can do that? I've looked at the integration tutorial, but can't find the heatmap in it.

See the DoHeatmap function in Seurat, which seems to be what that paper has used:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html#cb50
> DoHeatmap(pbmc, features = myGeneList, group.by="Condition")


text:
Splice variant analysis with ensembl genome/annotation/rna. Which files do I use?
I am running splice variant analysis. I wanted to use NCBI genome but the program works better with ensembl. I am a bit confused on primary vs. top level to use as my reference genome. I also do not know which annotation file to use. And I am not seeing a reference rna download. If anyone can help that would be immensely appreciated.

What is splice variant analysis?
For your reference question:
https://bioinformatics.stackexchange.com/questions/540/what-ensembl-genome-version-should-i-use-for-alignments-e-g-toplevel-fa-vs-p?rq=1
Generally using primary is the best option

Ensembl provides gene annotation files in both GTF and GFF3 format via the FTP site:
http://ftp.ensembl.org/pub/
The GTF and GFF3 files are named according to whether they contain data relating to the primary assembly, or on all top-level sequences including unplaced scaffolds, patches and haplotypes. More information can be found in the README files within each folder. Which genome sequence file and annotation file you use will depend on the details of your study and have been discussed in detail on other posts.
Ensembl also provides reference RNAseq data aligned to the reference genome sequence as BAM files with corresponding .bai index files on the FTP site:
http://ftp.ensembl.org/pub/current_bamcov/homo_sapiens/genebuild/
More information can be found in the README:
http://ftp.ensembl.org/pub/current_bamcov/homo_sapiens/genebuild/README


text:
Making a box plot of the proportion of cells in each cell types in two groups
I have number of cells in three cell types T, B, M in a Seurat object
For two groups of patients, cancers and controls like
> table(Idents(pbmc), pbmc$samples)

    PN0252_0001 PN0252_0002 PN0252_0003 PN0252_0004 PN0252_0005 PN0252_0006 PN0252_0008
  T        3550         568         136         268        1835           2        1678
  B         136         129          11          25         133           0         203
  M         120          22          18          26         159           0          59
> dput(table(Idents(pbmc), pbmc$samples))
structure(c(3550L, 136L, 120L, 568L, 129L, 22L, 136L, 11L, 18L,
268L, 25L, 26L, 1835L, 133L, 159L, 2L, 0L, 0L, 1678L, 203L, 59L
), .Dim = c(3L, 7L), .Dimnames = structure(list(c("T", "B", "M"
), c("PN0252_0001", "PN0252_0002", "PN0252_0003", "PN0252_0004",
"PN0252_0005", "PN0252_0006", "PN0252_0008")), .Names = c("",
"")), class = "table")
>

Patients PN0252_0001, PN0252_0002, PN0252_0003, PN0252_0004 are controls
and PN0252_0005, PN0252_0006, PN0252_0008 are cancers
I want to look at the difference of the proportion of cells in each cell type between cancer and control
But I am confused for getting proportion should I divide the number of cells on total cells in the object or total cells in each group
Basically I need a box plot like although here numbers are not correct

This looks like a ggplot with a facet added to it. What have you tried?

What I plotted here is what I have tried but I am not sure about the proportions

Your boxplot seems correct. What do you think is wrong about it?

I don't know proportion of T cells for example, I should divide the number of T cells in each patients on the total number of cells or total cells in controls and cancers separately

The answer to that depends on the question you're asking and what you wish to showcase. In any case, only you can answer it, we can't do that for you.

Ram RS as always your comments are extremely helpful . Thanks a million

@mahta-mira, are you able to write up your solution as an answer, then mark it as a selected answer? This will help the StackExchange bots work out what's going on.

At the moment I am not certain if my solution is correct :(


text:
Is it normal to have a smaller bam file after merging bam files?
I have 2 bam files that belong to the same sample and I merged them with samtools merge. And after merging, I realized that the merged version is a bit smaller than the sum of the other two separate files' size. I also wonder if I merged "the fastq" files of this sample first and convert it to bam would the result be the same?

Perhaps the input was in .sam format and output in .bam?

@user438383 No, I had "the fastq" files but before getting ".bam files" I didn't know they belong to the same sample. So I converted them to ".sam" and then to ".bam" respectively. And now that I learned that they actually belong to the same sample I have to merge them.

Converting fastqs to bams is technically doable, but why?  Do you mean merging the fastqs and realigning?

@swbarnes Yes, exactly.

If you merge a lot of BAM files you lose the overhead of the header which depending on the size of your BAM files this can be a significant difference or not.
With FASTQ files there should be less difference (as they don't have an overall header). However, the compression ratio can change depending on how you zip your files. [EDIT: I just tested it merging two fastq.gz that were 357MB and 562MB - the result was fastq.gz that was smaller than the sum at only 826MB]
No, the BAM header information would be lost in the conversion from BAM -> FASTQ -> BAM.


text:
Is it possible to add one more sample to the result of featureCounts?
I have some .bam files and I also got featureCounts by Subread. But after getting the featureCounts I realized that among my bam files 2 of them actually belong to the same sample. So I wonder if I can remove the result of these two and after merging bam files add them again into featureCounts without doing all the featureCounts for other samples again.

Sure, paste and company will help or load it into R and cbind stuff as needed, but seriously, featureCounts is blazingly fast and even with many samples on a modest laptop it just takes minutes. Just rerun it. Given the 59min period of time that this question is here it would have been finished by now probably ;-) For code suggestions it is mandatory to add example data though.

It sounds like you are supplying multiple BAMs to featureCounts, while this is possible I think it's better to run featureCounts on each BAM separately and then use a Python script or cbind function in R as mentioned above to combine the results of multiple samples. This shall give you more flexibility.

Yeah thanks for the answers. As ATpoint mentioned I realized that featureCounts doesn't need much computational power and time. And good to know I can merge with Python or R script.


text:
How to use multiple .pdb files as templates for homology modelling?
I just ran an AlphaFold job for PduK from the Pdu BMC operon and got very strange results including a large circular loop of random coil. According to Mayer et al., PduK should be semi-triangular protein that forms hexagonal homohexamers. To fix this issue, I am interested in performing homology modeling of PduK with other Pdu proteins (namely PduA, PduB, PduJ, PduK, PduU, and PduT, all of which have crystal structures in the Protein Data Bank) as templates. I was only able to figure out how to use one template .pdb file in I-TASSER with its "-restraint4" parameter and was not able to get MODELLER working for this purpose.
My PduK residue sequence is: VKQSLGLLEVSGLALAISCADVMAKAASITLVGLEKTNGSGWMVIKIIGDVASVQAAISTGVSFADQRDGLVAHKVISRPGDGILSHSVTPESESEPAPAPTPVVPHEEIPEDHAAPEAPQDAELISCNLCLDPACPRQKGEPRSLCLHSGKRGEA

Do you get the solution to this problem as I also want to do the modeling by using two templates at the same time?

AlphaFold2
You are asking about how to hybridise different threaded models, but I would not give up on AF just yet as it can do complexes.
AlphaFold2 is conceptually amazing, but does have some issues to such an extent that I blogged about common problems I was getting so many questions. Specifically, it models one protein chain in isolation.
In the RoseTTAFold paper David Baker and co. mimicked the AlphaFold2 algorithm hinted at from the CASP14 competition. They discovered that if multiple chains were given as a protein concatenation, with suitably long polyglycine linkers, a good oligomer or hetero complex could be resolved. Quickly after the public release of the Colabs notebooks of AlphaFold2 in the Twittersphere there were reports of homo/hetero-oligomers being resolved nicely. Of particular interest is the report that if you put a spurious pairing you get bad results —i.e. a control of sorts.
This strategy many be applicable to your homohexamer.
Your long random loop could be binding something else, so if you know what you may also then want to try with that something else!
Threading, ab initio and hybridisation
Given a template PDB and an aligned sequence a threaded model of the sequence can be made by replacing the residues where required. SwissModel, ModBase and Rosetta partial_thread app does this in the open, while Phyre2 and I-Tasser do it hidden in the pipeline.
If a span is absent, SwissModel and others will try to predict from scratch (ab initio) the span.
Features of different threaded models can be stuck together (hybridised) manually or automatically as do Phyre2 and I-Tasser. The RosettaCM (comparative modelling) toolkit is really handy in that one can customise how different parts are treated by adding Rosetta constraints (restraints in crystallography parlance), ligands etc etc. It should be noted that threaded models are faithful to the template, so can be converted into the same complexes as the template, but as soon as they get hybridised they loose this.
You currently have various threaded models against a panel of homologues of your choosing and are looking to hybridise them based on your requirements —I assume to circumvent a mistake the standard I-Tasser server makes, such as messing up oligomers.
So the hybridisation from RosettaCM (skipping the threading) is probably a good choice, however, requires familiarity with Rosetta.
The model would need to be "symmetric" (oligomer) before starting the hybridisation (Step 2 in the manual). But that can simply be a case of getting the oligomeric template and superimposing that many copies of your threaded model onto it. Then as your threaded model is from a different framework, it needs to be adapted to the forcefield in Rosetta which differs from I-Tasser, so needs to be minimised, using Relax.
Another minor note is that a weakness is the pairwise alignment for all of these, which for distant variants can be problematically erroneous —hence why Phyre2 and I-Tasser run a MSA and from there extract the pairwise, which is a lot more accurate. I do not know how the restraint4 parameter works.
Edit: SwissModel
I realised just now that maybe SwissModel was not tried —it is an old threader, but does oligomers nicely. In fact, I tried your sequence and the hexamer looks decent and the internal metric is pretty good...


text:
Sample versus subject in scRNA-seq
I am reading this paper, where the authors mention

We apply SAUCIE to the batch correcting, denoising and clustering of
an 11-million cell mass cytometry dataset with 180 samples from 40
subjects in a study of the dengue flavivirus

What exactly is the difference between a “sample” and a “subject” in scRNA-Seq?
PS: here’s a possibly relevant visual from the paper, but I still couldn’t understand the difference

A subject is an individual or donor. You can take multiple samples from the same subject, e.g. from different tissues, different timepoints etc.


text:
Beginner with DESeq2 having issue with analysis
I am trying to run DESeq2 analysis on my featureCounts output counts table.
I ran the following code:
library("DESeq2")

samplefile<-"C:/Users/RGB King/Documents/R/~/R/DESeq2_analysis/sample_table.txt"

coldata<-read.table(file=samplefile, sep="\t", header=FALSE, col.names = as.factor(c("","treatment")), row.names=1)

countfile<-"C:/Users/RGB King/Documents/R/~/R/DESeq2_analysis/feature_counts_test_AvABC.txt"
countdata<-read.table(file=countfile, sep="\t", header=FALSE, row.names = 1)

countdata<-countdata[-c(1:5)]

dds <- DESeqDataSetFromMatrix(
    countData = countdata,
    colData = coldata,
    design = ~ treatment
)

I then get the following error:
Warning message:

In DESeqDataSet(se, design = design, ignoreRank) :

  some variables in design formula are characters, converting to factors

The countdata table looks as follows:
|------|Variable1|Variable2|Variable3|Variable4|Variable5|Variable6|Variable7|Variable8|
|GeneID|A_1.bam--|A_2.bam--|A_3.bam--|A_4.bam--|AB_1.bam-|AB_2.bam-|AB_3.bam-|AB_4.bam-|
|THI3--|:12:-----|:34:-----|:8:------|:15:-----|:6:------|:22:-----|:41:-----|:33:-----|
|THI7--|:17:-----|:78:-----|:9:------|:23:-----|:11:-----|:54:-----|:19:-----|:7:------|

The coldata looks like this:
|--------|treatment|
|A_1.bam-|:A:------|
|A_2.bam-|:A:------|
|A_3.bam-|:A:------|
|A_4.bam-|:A:------|
|AB_1.bam|:B:------|
|AB_2.bam|:B:------|
|AB_3.bam|:B:------|
|AB_4.bam|:B:------|

So after I got the error I found posts on a forum that stated that DESeq2 only recognises "integer" class input values for the analyses. Therefore, I ran the following code:
sapply(countdata,class)

This then showed me that the data in the countdata table was in "character" class, so to fix this I ran the following:
chars <- sapply(countdata, is.character)
countdata[ , chars] <- as.data.frame(apply(countdata[ , chars], 2, as.integer))

But this then caused the values in the GeneID row to become "NA". So my questions are how can I change the class of just the GeneID row back to what it was? Will this then allow the analysis to run correctly? And in my code for the countdata I do specify "row.names =1" so why does this not keep my values in the GeneID row separate from the rest of the table?

If the first column is GeneID and your data.frame was called df then do df[,2:ncol(df)] or move GeneID to rownames and then convert to numeric class.

Hi there, @ATpoint thanks for taking the time to answer I appreciate it. Could you maybe explain how to move the GeneID to row.names in my code? If I can explain better, I need to just keep the cells in GeneID row i.e. "A_1.bam" as a "character" class and change the rest of the cells that contain the DE gene counts to "integer" class. I am not sure if this will work though, but I just want to try it and see if that is why I am getting the error I mentioned in post.

Cause when I run sapply(countdata,class) I still get the following output          V7          V8          V9         V10         V11         V12         V13         V14
"character" "character" "character" "character" "character" "character" "character" "character"

Why are you importing your count data with header = false when you clearly have a header?

Yes apologies I corrected that, however, I still can't resolve the issue of the data in the counts table I input being in "character" class. Is this a common error that one gets when working with featureCounts generated counts tables? It might be worth mentioning, I changed the file from what I got directly as the output file i.e. a file in ASCII text format to a text document in notepad on windows. I am not sure if that maybe caused DESeq2 to recognise the integers as "character" class as well? Would it be better to convert the output of featureCounts to a .csv file format in excel?

Are you sure you cant use the feature counts file without manipulating it in windows at all?  The warning you quote above (not an error) is about your coldata, not your count data.

@swbarnes2 I will give that a try then thank you. Could you maybe explain what is wrong with my coldata file please? And how should I correct it so that the analysis can run without the error?

There is nothing wrong with your colData. It is encoded as character rather than factor, and DESeq2 is turning it into factors, that is all that the warning (warning != error) tells. I never needed to convert my matrix to load featureCounts output into DESeq2. Just load the output with the header, use rownames() to move the gene column to rownames, remove all the non-relevant columns (that is usually 1-7) and then feed the remaining one into DESeq2. I have no featureCounts output at hand right now. If you need code then please provide suitable example data via the dput function for easy copy/pasting.

Hey @ATpoint sorry your suggestion on my other question reminded me about this. I got assistance from a friend who works as a bioinformaticist at a private company. He noted that there was an annotation step missing in my data and interestingly even with his code he wrote for me I still got the warning message I mentioned, which he attributed to using RStudio in Windows. I still would like to send examples of my data through to you though, as although I did get an output from his code I am unsure if it is correct.

As the warning message is suggesting, you should convert the class of your treatment column in coldata from character to factor.
Use this:
coldata$treatment = factor(coldata$treatment)

Thank you, I will give that a try


text:
shiny dashboard to visualise RNA seq data
I am trying to build a shiny dashboard for the visualization of different RNAseq data sets.
However, I am encountering problem with reactive input, the code is always using the second dataset that I have loaded/read in. Anyone has faced a similar issue? If yes, how did you resolve it.
My code is
    library(shinydashboard)
library(uwot)
library(DESeq2)
library(gridExtra)
library(tidyverse)
library(RColorBrewer)
library(DESeq2)
library(pheatmap)
library(DEGreport)
library(vsn)
library(RColorBrewer)
library("genefilter")
library(org.Hs.eg.db)
library(dplyr)
library(tidyverse)
library(fgsea)
library(clusterProfiler)
library(ggplot2)

set_1<-load("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData2.RData")
set_2<-load("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData1.RData")
data_list = list(set_1=set_1,set_2=set_2)

ui <- dashboardPage(
    dashboardHeader(title = "Data Visualizer", titleWidth = 300),
    dashboardSidebar(
        width = 300,
        sidebarMenu(

            menuItem("Datasets",  icon = icon("cog"),
                     selectInput("Datasets", "Datasets:", choices = list("sample1" = "set_1", "sample2" = "set_2"),
                                 selected = "set_1")),
            menuItem("Quality Control", tabName = "widgets", icon = icon("th")),
            menuItem("Differential Genes", tabName = "widgets2", icon = icon("th")),
            menuItem("Downstream", tabName = "widgets3", icon = icon("th"))

         )
    ),
    ## Body content
    dashboardBody(
        tabItems(
            # First tab content
            # Second tab content
            tabItem(tabName="widgets",
                   h2("widgets"),
                   plotOutput("widgets"),

            ),
            tabItem(tabName = "widgets2",
                    h2("Widgets2 tab content"),

            ),
            tabItem(tabName = "widgets3",
                    h2("Widgets3 tab content"),
                    plotOutput("widgets3"),

            )

        )
    )
)

server <- function(input, output) {

    datasetInput <- reactive({
        df <- data_list[[input$Datasets]]
        })

        output$widgets <- renderPlot({
            datasetInput()
            par(mfrow=c(1,2))

            boxplot(counts(dds, normalized=F), outline=F, col=dds$condition, medcol = "white", cex.axis=0.6, main="Without Normalization")
            boxplot(counts(dds, normalized=T), outline=F, col=dds$condition, medcol = "white", cex.axis=0.6, main="Normalized")

            })

        }

shinyApp(ui, server)

There are two issues that I can see.
The first is that the results are loaded via RData files, which include information about the variable names. Assuming these work, I expect they create a variable within the R environment called dds which has the DESeq2 results. The second load overwrites that variable:
set_1<-load("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData2.RData")
set_2<-load("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData1.RData")

What should probably be done is storing the results as an R Dataset that doesn't include variable name information:
# create set_1 DESeq2 results called 'dds'
saveRDS(dds, "C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData2.rds")
# create set_2 DESeq2 results called 'dds'
saveRDS(dds, "C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData1.rds")

Then those can be loaded into the R environment using the assignment you expect:
set_1<-readRDS("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData2.RData")
set_2<-readRDS("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData1.RData")
data_list = list(set_1=set_1,set_2=set_2)

The second issue is that the plot output displays dds, rather than the datasetInput variable, and there is no assignment to change the definition of dds:
    output$widgets <- renderPlot({
        dds <- datasetInput() ## Note: added assignment
        par(mfrow=c(1,2))

        boxplot(counts(dds, normalized=F), outline=F, col=dds$condition, medcol = "white", cex.axis=0.6, main="Without Normalization")
        boxplot(counts(dds, normalized=T), outline=F, col=dds$condition, medcol = "white", cex.axis=0.6, main="Normalized")

        })

    }

Thank you, your tips were very helpful


text:
How to know if the DNA sequence has been assembled and why is it important to know how it was assembled?
I have downloaded my FASTA format files, that have the DNA sequences of the coding region of the genes and the DNA sequence of the complete genome, from NCBI. How can I recognize if these sequences have been assembled or not ?
In particular, I read this paper Genome Sequencing of a Genetically Tractable Pyrococcus furiosus Strain Reveals a Highly Dynamic Genome and I downloaded the sequences that they indicate in Nucleotide sequence accession number where they say:

The fully assembled DNA sequence of COM1 has been deposited in GenBank (accession number CP003685).

I clicked on the link CP003685 and I downloaded the DNA sequences that I mentioned at the beginning by clicking on "Send to", then "Complete Record",then "File" (in Choose Destination), then FASTA (in Format). For the cds file I made the same steps but clicking on "Coding Sequences" and I could do the same to have the sequences of the genes clicking on "Gene Features". But are these sequences assembled? By reading the extract of the article I would say yes... but is there a way to know this from the FASTA files that I downloaded or the NCBI website ?
Moreover, why is it important to know how a sequence has been assembled (so which software has been used) ?
I apologize for my very basic questions and I thank you in advance.

Genomes are always assemblies. Assembly basically means making a consensus sequence from many partially overlapping readouts. As there is currently no sequencing technology that allows to sequence an entire genome or even an entire chromosome with a single super high fidelity read that is without errors we have to combine several reads to form the consensus which is then the assembly.

Moreover, why is it important to know how a sequence has been assembled (so which software has been used) ?

Tbh, I never cared on how a genome I used was assembled. To be fair I only work with mouse or human, and for people who work with genomes that are not as well curated as these it might indeed matter, but others would need to comment on that, not really my field.

Thank you @ATpoint ... so by "consensus" you mean that you find the nucleotide that appears most frequently at each position by looking at the overlapping readouts? Thus when you have a fasta format file that starts, for example, with ">NC_003413.1 Pyrococcus furiosus DSM 3638, complete sequence" and then you have the nucleotides sequence, does it always represent an assemblied genome ?

I guess so, yes.

A nucleotide sequence may be a part of a genome assembly. If so, just click the assembly link from the right panel and you'll be redirected to the assembly page.
In this assembly page, it says the assembly only contains one chromosome.
Assembly method is also available for some assemblies.


text:
Question about the practical relevance of chromosome position and p-value
Article of interest: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3023908/#SD1
Referring to Figure 2b of that article, if I understand correctly, the x-axis refers to the position of each variant in a chromosome while the y-axis is a log transform of the raw p-value of association of each variant with some variable of interest.
My questions is, is there practical importance in finding strong association (high log transform p-value) around one particular position versus finding it on some other position in the same chromosome?
My follow up question is, are these log transformed p-values publicly available?
Context: I am a statistics student studying multiple testing. My biology background is very limited.

In the specific case that you are talking about, the answer is very much yes that the specific position is important. Genes occur at (fairly) fixed positions along chromosomes.
The point of that figure, as indicated in the legend, is that the very low p-values are around the gene RPM1. That gene is known to have important natural genetic variation associated with the phenotype in question (the one for whichg they are looking for genetic associations). That figure is therefore an important validation of the methods of the authors.
In terms of finding such data, a fun account on twitter is the GWASbot, that just posts random manhattan plots (which is the informal name for the plot you are talking about), along with links to source data. I believe that that source data includes the p-values for each ascertained position, along with other summary statistics. One example is here.
If you are specifically interested in the topic of multiple test correction, an important related concept in GWAS is that of genome-wide significance.
Update:
To unzip and look at the file, use bgzip:
$ bgzip -d ~/Downloads/continuous-20491-both_sexes.tsv.bgz
$ head ~/Downloads/continuous-20491-both_sexes.tsv
# prints the following:
chr pos ref alt af_meta beta_meta   se_meta pval_meta   pval_heterogeneity  af_AFR  af_CSA  af_EUR  beta_AFR    beta_CSA    beta_EUR    se_AFR  se_CSA  se_EUR  pval_AFR    pval_CSA    pval_EUR    low_confidence_AFR  low_confidence_CSA  low_confidence_EUR
1   11063   T   G   2.088e-02   2.247e-02   2.072e-01   9.137e-01   0.000e+00   2.088e-02   NA  5.212e-05   2.247e-02   NA  1.287e-01   2.072e-01   NA  5.080e-01   9.137e-01   NA  8.000e-01   false   NA  true
1   13259   G   A   2.805e-04   9.840e-02   8.909e-02   2.694e-01   1.000e+00   NA  NA  2.805e-04   NA  NA  9.840e-02   NA  NA  8.909e-02   NA  NA  2.694e-01   NA  NA  false
1   17641   G   A   8.908e-04   -1.244e-02  4.922e-02   8.005e-01   1.000e+00   NA  NA  8.908e-04   NA  NA  -1.244e-02  NA  NA  4.922e-02   NA  NA  8.005e-01   NA  NA  false
1   30741   C   A   NA  NA  NA  NA  NA  9.400e-04   NA  NA  3.935e-03   NA  NA  1.297e+00   NA  NA  9.976e-01   NA  NA  true    NA  NA
1   51427   T   G   NA  NA  NA  NA  NA  1.127e-03   NA  NA  8.780e-01   NA  NA  9.366e-01   NA  NA  3.486e-01   NA  NA  true    NA  NA
1   57222   T   C   7.078e-04   -4.153e-03  5.656e-02   9.415e-01   1.000e+00   NA  NA  7.078e-04   NA  NA  -4.153e-03  NA  NA  5.656e-02   NA  NA  9.415e-01   NA  NA  false
1   58396   T   C   2.735e-04   3.011e-02   8.221e-02   7.142e-01   1.000e+00   2.622e-03   NA  2.735e-04   -1.069e-01  NA  3.011e-02   7.079e-01   NA  8.221e-02   8.800e-01   NA  7.142e-01   true    NA  false
1   62745   C   G   NA  NA  NA  NA  NA  2.701e-03   NA  NA  5.524e-01   NA  NA  5.529e-01   NA  NA  3.178e-01   NA  NA  true    NA  NA
1   63668   G   A   NA  NA  NA  NA  NA  NA  1.129e-03   3.673e-05   NA  2.205e-01   -4.462e-02  NA  9.738e-01   2.544e-01   NA  8.209e-01   8.608e-01   NA  true    true

Thank you. I tried downloading an item from the GWAS bot. It was a .bgz file. May I ask what the best way is to extract data from this filetype? I tried to open it with Excel but it did not work.

@variancekills See update-  you need to decompress the file, then it will be a TSV that you can in theory look at in excel, though it is pretty big for excel (several GB). I would suggest R, python or similar if you are familiar

R is fine for me. Thank you very much for your help! I've been working on a small pet project for multiple testing that requires a very specific data setup. I need thousands of hypotheses that I can line up meaningfully. This setup works very well for testing my method.

thank you very much. i was able to open the file. May I ask you one more question? Which column in the example you posted is the actual p-value used to compute the logtransformed p-values on the plot?

Several of these columns represent p-values. All of the column headers beginning "pval" are such, for different ancestry groups (e.g. pval_EUR). For your purposes, I believe that the overall p-value (pval_meta) is probably sufficient.

Thanks. Yes, I was able to reconstruct the map using pval_meta. This is very helpful. Thank you very much!

Are these log transformed p-values publicly available?

For human population data, it's pretty common for the association statistics to be either included with the research paper, or uploaded to a public database (e.g. see curated compilations like GWAS Central or NHGRI GWAS Catalog).
For plant and other non-human data, there doesn't seem to have been as much drive to include GWAS data with the research paper, and it's often up to the authors to decide to do that. I find this a little odd, given that there should be fewer ethical barriers to releasing non-human public data.
If you're interested in a particular dataset, you can contact the authors, but I acknowledge that such an approach can be too high of a barrier for most interested researchers.

Is there practical importance in finding strong association (high
log transform p-value) around one particular position versus
finding it on some other position in the same chromosome?

It's common for statistical tests for association to assume independence of variant locations. With this in mind, two associated regions that are one base apart should not be any more interesting than two associated regions that are 50 megabases apart.
The problem with this assumption is that DNA variation is not statistically independent. Recombination of DNA (and associated reassortment of genetic variation) between two locations is more common when the locations are close, and this linkage should really be accounted for when doing association tests. One way to do this is to combine variants that share the same ancestral origin into a haplotype, then do association tests on the haplotype block rather than the single-base genotype. Even if the particular observed variants are not directly causative, there may be another untyped causative variant (or combination of variants) on the same haplotype block, meaning that the block becomes a proxy for the causative variant.
But this haplotype-based approach becomes quickly computationally intractable when looking at larger regions of chromosomes, and is very difficult to use for comparing populations with substantially different ancestral backgrounds. The workaround (if it could be called that) is to ignore the linkage aspect of genetic variation, and return to those same statistical tests, but declare regions that have peaks of close association to be statistically interesting, and zoom into them to look in more detail at the region.
That zoom in would be a great opportunity to use more complex / accurate haplotype-based tests for confirming an associative signal around a particular region. But for the most part, as with the paper you've indicated, the tests are identical.
However, even ignoring the haplotype issue... <begin rant>
The most common manhattan plots represent only p-value on the Y-axis. I think this is a mistake, because the p-value only represents the confidence in a statistic, and the important thing for association is the statistic itself, rather than the confidence in that statistic.
I have a few specific concerns about this pure p-value approach:

there may be values that have less confidence, but a much higher associative statistic, and contribute more to the trait of interest (e.g. see here)
values with a low association statistic are more sensitive to genotyping error, leading to spurious associations.
the p-value doesn't indicate the direction of association
p-values should not be used for ranking, and this ranking is implicitly done when using the p-value only manhattan plot

In addition, the generated p-values often don't make sense, and I haven't seen any commentary on this. Anything smaller than $10^{-40}$ doesn't represent any reasonable natural random occurrence, even on galactic timescales, but often these astronomically-small p-values can become insignificant through a few random mutations in one of the populations.
As a quick-fix using the existing statistics that are typically presented in GWAS results, where single points are desired, I recommend that people instead use the standard deviation to determine the minimum expected beta score:

see here for more details.
More rigorously, I think that people should consider doing their own internal validation of associations via bootstrap sub-sampling, and report a non-parametric statistic: the best observed rank after multiple bootstrap sub-samples. This requires access to genotype-level data, so cannot be done on existing GWAS results. I explore this in more detail in this preprint.
Yes, I'm aware this goes against the grain of how GWAS is typically done. I have been emboldened in my position on this after reading the ASA principles about the p-value, and seeing that many of the principles are in direct opposition to the typical manhattan plot:
https://www.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#_i31

I think that this is a fair, well-argued criticism of GWAS methodology, supported by references, and which I personally tend to agree with. I myself have serious problems with NHST, specifically in high-n cases like genomics. But I'm not sure that it really gets at the asker's question regarding getting a baseline understanding of what GWAS is and some of the fundamental ideas that motivate it. And yet, I do think that it's valuable to have this answer's viewpoint represented on this page.

Yes, fair comment. I've added a couple of direct answers to the specific questions asked. My additional rant was on a related issue which answers other questions which I think are more relevant for interpreting GWAS data.


text:
Getting all promoters of a specific gene
I'd like to perform a motif analysis of all promoter regions of the gene ENSMUSG00000020538 in mus musculus.
To do so, I wanted to use Biomart :
library('biomaRt')
mart <- useMart("ensembl")
mart <- useDataset("mmusculus_gene_ensembl", mart = mart)
all.genes <- getSequence(id = 'ENSMUSG00000020538',
                         type="ensembl_gene_id",
                         seqType="coding_gene_flank",
                         upstream=100,
                         mart=mart)

But this provide me only the flanking region of the sepcified gene, not all its promoter regions. How could I do ?

Answer from @ATpoint, converted from comment:
You need to extract the TSS of all transcripts of that gene, not the gene borders. Query biomaRt for that gene, then set options to extract all transcripts, then get the sequence for the TSS coordinates of those (with a window).


text:
Counting the co-occurrence of variants A and B in an aligned sequencing read
I need some help getting started on this project.
To simplify we want to be able to quantify the occurrence of 3 variants on each sequencing read in an alignment file as a proxy measurement for recombination.
The possible combinations we have are A, B or AB.
For example, variant A is at position 100 and variant B is at position 300 of the reference and AB would have both.
We have alignment files against a reference sequence, and then I get stuck. I could technically look by eye in tablet or IGV for example but is there a way to extract each read then query by expected base at the reference position then bin these reads into the categories; A, B, AB and unclassified then calculate the total occurrences of each?
Thanks for your time and suggestions!

You'll need to write something to do this, using pysam or a similar package.

Related question: Access base aligned to particular reference position
After an initial effort posted here, I ended up rewriting this and testing it a little, and posted the script here. It still has some drawbacks (handling indels intelligently) but now it more or less works.


text:
FoldX PositionScan functions only works for certain position inputs
I'm using the PositionScan function in FoldX 5.0 to calculate delta delta G for some mutants. I'm working with the command prompt in Windows 10. For some reason, this function will only return output for residues 1-9 of my protein - any residues in the double and triple digits don't work.
For example, my protein only has one chain (A) and it has Leu at position 9 and Lys at position 10. For whatever reason, the following command will work:
FoldX --command=PositionScan --pdb=protein.pdb --positions=LA9a

While this one won't:
FoldX --command=PositionScan --pdb=protein.pdb --positions=KA10a

I have triple checked the PDB file itself to make sure I'm not typing in the wrong letter codes and positions. The format of the PDB file seems to be fine, when I open it in a structure viewer nothing seems out of the ordinary.
What's going on here?
Thanks for any help in advance! I'm new at this stuff.

What happens when you try other commands?
I see and appreciate that you triple checked the PDB file, but it totally sounds like a problem there —e.g. the PDB file is tabulated incorrectly (cf. PDB ATOM definition, hence the digit issue.
Say what does the PDB look like after Optimize command (energy minimisation) or PDBFile (a simple roundtrip)? If you get a nonapeptide, then you need to quadruple check!


text:
TCR-seq or scRNA-seq
I need idea, intuition, suggestion please
I have 10X scRNA-seq of PBMC (multiome's 3' poly-A capture), can I capture TCR/BCR from my scRNA-seq or I need TCRseq ? I googled and likely TCR and BCR immune profiling require 5' transcript capture and the V(D)J variable sequence is on the 5' end of the gene. So the 3' capture typically yields the C constant region of TCR and BCR transcripts. can you help me please? Thank you

The 3' gene expression protocol will capture TCR and BCR mRNAs but this may not be very helpful to you. As you already mentioned only the 3' end will be sequenced, which are the constant regions. With this data you may potentially be able to discriminate TRAC from TRBC1 and so forth.
However, in most cases researchers are interested in clonotypes (combinations of variable regions) and these you can only obtain by using the dedicated TCR or BCR kits from 10X. They contain specific sets of primers for these genes to cover the whole TCR or BCR gene body and not only the 3' end. This way you can see the different clonotypes in your population
The human TCR loci can be found for example at HUGO. The last gene at each locus is the constant and is identified by a C at the end. TRAC = T cell receptor alpha constant (I did not use the correct naming convention in the original answer, sorry)
B-cell receptors basically are membrane bound immunoglobins (IG). Therefore you will likely find the constant regions of the heavy chain in your data. The different heavy chains are what identifies the IG. For example IGHA1 is the heavy chain of an alpha IG. If the secreted and membrane bound form of the IGs is always the same per B-cell, unfortunately, is beyond my knowledge of B-cells.

Thank you I got a question, in this data base I could not find anything corresponding to TCRVδ2 , you think in which part I can find this gene?

The nomenclature used by immunologists may not match the official gene names (which are in HUGO). The TCR gamma variable genes are named TRGV followed by a number. Problems with different names for the same gene / protein are all over the place and you will sometimes have to search several data bases (HUGO, UniProt, Entrez ... ) or the literature for an answer

According to our experience 3’ sequencing does not produce enough data to recover TCR and BCR repertoire. But we have some experience recovering repertoires from non enriched 5’ 10x RNA-Seq GEX data.
With MiXCR we have specific presets that allow analyse 10x data in one single line command:
mixcr analyze 10x-5gex-full-length \
 --species hsa \
 sample_R1.fastq.gz \
 sample_R2.fastq.gz \
 sample_result

Also feel free to check out our documentation portal:
https://docs.milaboratories.com/mixcr/reference/overview-built-in-presets/#10x


text:
loading in .txt files with different spacing in column names into R
Trying to read in this file but saved as a .txt file: http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt
gene_conversion = read.table("HOM_MouseHumanSequence.rpt.txt",header=T, fill=T)
but the column names are splitting weirdly.
> head(gene_conversion)
        DB  Class        Key Common Organism       Name        NCBI Taxon    ID                      Symbol
1 38710983 mouse, laboratory  10090     Gdnf      14573  MGI:107430 Chr15   3.8                          cM
2 38710983  human       9606   GDNF     2668  HGNC:4232 OMIM:600837  Chr5 p13.2   Chr5:37812677-37840044(-)
3 38710984 mouse, laboratory  10090     Rrm1      20133   MGI:98180  Chr7 54.72                          cM
4 38710984  human       9606   RRM1     6240 HGNC:10451 OMIM:180410 Chr11 p15.4    Chr11:4094685-4138925(+)
5 38710985 mouse, laboratory  10090   Hoxd12      15432   MGI:96204  Chr2 44.13                          cM
6 38710985  human       9606 HOXD12     3238  HGNC:5135 OMIM:142988  Chr2 q31.1 Chr2:176099802-176100760(+)
                                                            EntrezGene
1                                             Chr15:7840327-7867056(+)
2 NM_199231,NM_001278098,NM_199234,NM_001190469,NM_001190468,NM_000514
3                                          Chr7:102090902-102118978(+)
4                     NM_001033,NM_001318065,NM_001330193,NM_001318064
5                                            Chr2:74505357-74508049(+)
6                                                            NM_021193
                                                                                ID.1
1                                     NM_001301357,NM_010275,NM_001301333,NM_001301332
2 NP_001265027,NP_000505,NP_001177397,NP_001177398,XP_016864826,XP_011512332,NP_954701
3                                                                            NM_009103
4                                     NP_001024,NP_001317122,NP_001304994,NP_001304993
5                                                                            NM_008274
6                                                                            NP_067016
                                             Mouse    MGI ID.2 HGNC ID.3 OMIM Gene ID.4 Genetic Location Genomic
1 NP_001288286,NP_001288262,NP_034405,NP_001288261 P48540   NA   NA   NA   NA   NA   NA      NA       NA      NA
2                                           P39905          NA   NA   NA   NA   NA   NA      NA       NA      NA
3                                        NP_033129 P07742   NA   NA   NA   NA   NA   NA      NA       NA      NA
4                                           P23921          NA   NA   NA   NA   NA   NA      NA       NA      NA
5                                        NP_032300 P23812   NA   NA   NA   NA   NA   NA      NA       NA      NA
6                                           P35452          NA   NA   NA   NA   NA   NA      NA       NA      NA
  Coordinates X.mouse. X. human. X..1 Nucleotide RefSeq IDs Protein RefSeq.1 IDs.1 SWISS_PROT IDs.2
1          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
2          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
3          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
4          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
5          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
6          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA

I see spaces in the headings, which suggests this might be a tab-separated file, rather than a space-separated file. If that's the case, you should be using sep="\t" in your read.table command, or alternatively using read.delim:
gene_conversion <-
  read.table("HOM_MouseHumanSequence.rpt.txt", sep="\t", header=T, fill=T)

or:
gene_conversion <- read.delim("HOM_MouseHumanSequence.rpt.txt")

On a related note, how can I subset based on just a column?

Can you please ask a separate question for that?


text:
converting mouse genes to human genes
I'm trying to convert mouse genes from PTX data to human genes in order to do a comparison with patient data to see what genes are being conserved. I'm using this file for the orthologs.
http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt
I read the table with:
gene_conversion <- read.table("HOM_MouseHumanSequence.rpt.txt", sep="\t", header=T, fill=T)

I have subsetted the mouse genes from the PTX seurat object to only include the mouse genes with:
mouse_only <- subset(scData, features = mouse_1)
head(rownames(mouse_only))
[1] "Xkr4"    "Gm37180" "Gm37363" "Gm19938" "Sox17"   "Gm6085"

Now I'm not sure how to convert the mouse genes in my mouse-only seurat object to match the human genes in the file.
Any help would be greatly appreciated.

Just checking / confirming... shouldn't that be a double equals in the subset command? are you using features == mouse_1?

No, just a single equals. When I run the command with double I get an error.
subset(scData, features == mouse_1) Error in FetchData(object = object, vars = unique(x = expr.char[vars.use]),  :    None of the requested variables were found:

Oh, I see now! Thanks. This seems to be an undocumented / new feature of subset that I haven't encountered before.

Any suggestions on how to continue?

Mice genes are usually the same as the human ones, but lower case with only the first character as upper (Sox17 -> SOX17).
toupper(rownames(mouse_only)) will do the trick.
Unfortunately, there are exceptions, and it is always better to relay on information from databases. The function bellow takes a vector of human genes, and uses biomaRt to return the corresponding genes in mouse.
# from Human to Mouse
convertHumanGeneList <- function(x){

   library("biomaRt")
   human <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
   mouse <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

   genesV2 <- getLDS(attributes = c("hgnc_symbol"), filters = "hgnc_symbol",
                 values = x , mart = human, attributesL = c("mgi_symbol"), martL = mouse, uniqueRows=T)

   humanx <- unique(genesV2[, 2])

   return(humanx)
}

You simply need to invert the 2 and you should be able to convert your vector rownames(mouse_only)
#---- Update using SOX17 as example

That didn't work. genesV2 is empty.

Works for me. I added a figure in my answer to show you the output I get

Figured out the issue, I didn't assign x a value. Now that I have the human genes, how can I do a correlation between the mouse genes from the PTX data vs the genes from patient data?

That's probably a topic for another question :D. You should accept the answer (if you managed to convert the genes name), and open a different topic.

library(dplyr)
#create the function to convert mouse to human genes
homologs = read.csv("http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt",sep="\t")
ms2hs <- function(gene_list){
  output = c()
  for(gene in gene_list){
    class_key = (homologs %>% filter(Symbol == gene & Common.Organism.Name=="mouse, laboratory"))[['DB.Class.Key']]
    if(!identical(class_key, integer(0)) ){
      human_genes = (homologs %>% filter(DB.Class.Key == class_key & Common.Organism.Name=="human"))[,"Symbol"]
      for(human_gene in human_genes){
        output = append(output,human_gene)
      }
    }
  }
  return (output)
}

#apply the function
msgenes=c("Zfp192", "Ifitm2", "Mettl8", "Kcnh2", "Npdc1", "Pitpnc1", "Ooep", "Mettl8")
ms2hs(msgenes)


text:
How to use SRA Toolkit to blast SRA accessions?
I'm trying to look for a gene in a group of SRA files. Web BLAST can't cope with them, so I assume they're too large. I've tried tblastn_vdb from the SRA toolkit, as follows. It had previously worked well for me for WGS files.
module load anaconda3/personal
source activate sratoolkit

cd ~/sra-tools
PATH=$PATH:~/sra-tools/sratoolkit.2.10.5-ubuntu64/bin
tblastn_vdb -db "SRR6256958" -query ~/sra-tools/query/query.txt -out ~/sra-tools/out/results.txt

I get nothing in my output file except for 'TBLASTN 2.6.1+'. My error file goes as follows:
Loading anaconda3/personal
  Loading requirement: fix_unwritable_tmp fix_setxattr
VDB: 2021-07-16T12:03:22 . debug: KSysDirResolvePath_v1 = '/[path]/sra-tools/sratoolkit.2.10.5-ubuntu64/schema'
VDB: 2021-07-16T12:03:22 . debug: KSysDirResolvePath_v1 = '/tmp/pbs.4070476.pbs'
VDB: 2021-07-16T12:03:22 . debug: KSysDirResolvePath_v1 = '/tmp/pbs.4070476.pbs'
BLAST engine error: NCBI C++ Exception:
    T0 "/panfs/pan1/sra-test/vdb3/git/blast_vdb/sratoolkit_2_9_0/c++/src/algo/blast/api/prelim_stage.cpp", line 304: Error: BLAST::ncbi::blast::CBlastPrelimSearch::Run() - search cannot proceed due to errors retrieving sequences from databases

/var/spool/PBS/mom_priv/jobs/4070476.pbs.SC: line 7: 4032551 Aborted [then it lists the original code again]

Does anyone know what I am doing wrong? Help would be much appreciated!

When you say "SRA files", are these fastq? If these files are 'large', BLAST may not be the right tool for this

@Chris_Rands No particular format at the moment since I haven't downloaded them separately. It was supposed to happen automatically? There should be enough space at least for this first file (8.79 gigabases), but the download doesn't even start.

SRA contains raw sequencing reads, blast usually serves in aligning genes / parts of genomes / assembled sequences (not sequencing reads). But I don't know tblastn_vdb. You sure the command is supposed to do what you think it should?


text:
calculating kmer nucleotide frequency per column
I have a list of sequences:
CAGGTAGCC
CCGGTCAGA
AGGGTTTGA
TTGGTGAGG
CAAGTATGA
ACTGTATGC
CTGGTAACC

Each sequence is nine nucleotides long. I want to calculate the dinucleotide frequency of AA, CG ..etc for each of the dinucleotide positions.
Similar problem shared with one nucleotide in this post.
For example, 1st two column will have these bases CA,CC,AG,TT etc, and second two should be AG,CG,GG,TG and so on. Output should be something like this
AA   0.25    0.34    0.56    0.43    0.00    0.90    0.45    0.34    0.31
CC   0.45    0.40    0.90    0.00    0.40    0.90    0.30    0.25    0.2
GC   0.00    0.00    0.34    1.00    0.30    0.30    0.35    0.90    0.1
TC   0.24    0.56    0.00    0.00    1.00    0.34    0.45    0.35    0.36

AA, CC, GC and TC could be any dinucleotide combination.
What I tried from the last post: mainly taken from @Devon
import itertools
hl = []

#making all possible dinucleotide combinations
bases=['A','T','G','C']
k = 2

kmers = [''.join(p) for p in itertools.product(bases, repeat=k)]
kmers = dict.fromkeys(kmers, 0)
print(kmers)
for i in range(9):
    hl.append(kmers)

# my seq file
f = open("H3K27ac.seq")
nLines = 0

for line in f:
    id = 0
# trying to read dinucleotide at each loop
    for (fst, sec) in zip(line[0::2], line[1::2]):
        id = id + 1
        hl[id][fst+sec] += 1
    nLines += 1
f.close()

nLines = float(nLines)
for char in ['AA', 'TT', 'CC', 'GG', 'CG', 'GC']:
    print("{}\t{}".format(char, "\t".join(["{:0.2f}".format(x[char]/nLines) for x in hl])))

By somehow it's estimation the total count of all dinucleotide in whole file.
I am trying to make it more general so that it can be used for tri- and any further number of nucleotides.

This should do the trick. The following code works for any aligned character set, but assumes that all sequences should be the same length.
import itertools

with open("f.txt") as seqs:
    sequences = list(line.strip() for line in seqs.readlines())

# Just in case sequence lengths are different
n = min(len(s) for s in sequences)
if any(n != len(s) for s in sequences):
    print("Sequence lengths are different!")

# Get all bases.
bases = set(c for seq in sequences for c in seq)

k = 2
kmers = [''.join(p) for p in itertools.product(bases, repeat=k)]
pos_to_freq = []

for i in range(n - k + 1):
    kmers_dict = dict.fromkeys(kmers, 0)
    for s in sequences:
        kmers_dict[s[i:i+k]] += 1/len(sequences)
    pos_to_freq.append(kmers_dict)

for kmer in kmers:
    print("\t".join([kmer] + [str(round(d[kmer], 3)) for d in pos_to_freq]))

Output for your example is below (I apologize for the formatting):
AA  0   0.143   0   0   0   0.143   0   0
AT  0   0   0   0   0   0.286   0   0
AC  0.143   0   0   0   0   0   0.143   0
AG  0.143   0.143   0.143   0   0   0.143   0.286   0
TA  0   0   0   0   0.571   0   0   0
TT  0.143   0   0   0   0.143   0.143   0   0
TC  0   0   0   0   0.143   0   0   0
TG  0   0.286   0.143   0   0.143   0   0.429   0
CA  0.286   0   0   0   0   0.143   0   0
CT  0.143   0.143   0   0   0   0   0   0
CC  0.143   0   0   0   0   0   0   0.286
CG  0   0.143   0   0   0   0   0   0
GA  0   0   0   0   0   0.143   0   0.429
GT  0   0   0   1.0 0   0   0   0
GC  0   0   0   0   0   0   0.143   0.143
GG  0   0.143   0.714   0   0   0   0   0.143

Works perfect with python3. Many thanks