text:
Run gmt format download from GSEA
I built the matrix as the gmt example in https://docs.gsea-msigdb.org/#GSEA/Data_Formats/ but It still has errors when I try to run it.
When I try the gene set from https://www.gsea-msigdb.org/gsea/msigdb/human/geneset/HALLMARK_ADIPOGENESIS.html
in .gmt format and load it with GSEA appears this error:
Thank you so much
Nerea

---- Full Error Message ----
col:34 > matrix's fColCnt:34

---- Stack Trace ----
# of exceptions: 1
------col:34 > matrix's fColCnt:34------
java.lang.ArrayIndexOutOfBoundsException: col:34 > matrix's fColCnt:34
at org.gsea_msigdb.gsea/edu.mit.broad.genome.math.Matrix.getColumnV(Matrix.java:261)
at org.gsea_msigdb.gsea/edu.mit.broad.genome.objects.DefaultDataset.getColumn(DefaultDataset.java:289)
at org.gsea_msigdb.gsea/edu.mit.broad.genome.objects.TemplateFactory.extract(TemplateFactory.java:85)
at org.gsea_msigdb.gsea/edu.mit.broad.genome.alg.DatasetGenerators.extract(DatasetGenerators.java:292)
at org.gsea_msigdb.gsea/edu.mit.broad.genome.alg.DatasetGenerators.extract(DatasetGenerators.java:288)
at org.gsea_msigdb.gsea/xtools.gsea.AbstractGsea2Tool.execute_one(AbstractGsea2Tool.java:79)
at org.gsea_msigdb.gsea/xtools.gsea.AbstractGsea2Tool.execute_one_with_reporting(AbstractGsea2Tool.java:103)
at org.gsea_msigdb.gsea/xtools.gsea.Gsea.execute(Gsea.java:166)
at org.gsea_msigdb.gsea/edu.mit.broad.xbench.tui.TaskManager$ToolRunnable.run(TaskManager.java:391)
at java.base/java.lang.Thread.run(Unknown Source)

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

Please show the first few lines of your input file; from the error, this looks like it might be a formatting issue, and without more information about the format of your file we can't help out.

Hello It was a formatting issue. Sorry, I solved the problem. Thank you

Great! glad it's resolved. Could you kindly add your resolution as an answer please? Firstly it helps the site because unresolved question continue to circulate, secondly it helps your reputation on the site. This is useful should you wish to ask another question here in the future.


text:
How can I change all text in a .tsv file to 1, except the first three fields?
I want to convert a .tsv file table in binary format (1 and 0).  In which all cells containing text would become 1 and blank cells would become 0,  except first three columns.
The given .tsv file has 6 fields
head adrbvbeT_copy.tsv
chr1    1   100
chr1    101 200
chr1    201 300
chr1    301 400
chr1    401 500
chr1    501 600
chr1    601 700
chr1    701 800
chr1    801 900
chr1    901 1000

tail adrbvbeT.tsv
chrY    57226501    57226600
chrY    57226601    57226700    Mmusculus-jaspar2016-Tcf21-MA0832.1
chrY    57226701    57226800    Mmusculus-HOCOMOCOv10-TF7L2_MOUSE.H10MO.A
chrY    57226801    57226900    Mmusculus-HOCOMOCOv10-STAT1_MOUSE.H10MO.A
chrY    57226901    57227000    Mmusculus-HOCOMOCOv10-ZEP2_MOUSE.H10MO.D
chrY    57227001    57227100    NA-SwissRegulon-SMARCC2.SwissRegulon    Mmusculus-HOCOMOCOv10-PBX2_MOUSE.H10MO.C
chrY    57227101    57227200    Mmusculus-HOCOMOCOv10-ZIC3_MOUSE.H10MO.C    Mmusculus-HOCOMOCOv10-TBX5_MOUSE.H10MO.D    Mmusculus-HOCOMOCOv10-PIT1_MOUSE.H10MO.C
chrY    57227201    57227300    Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D    Mmusculus-UniPROBE-Tbp.UP00029  Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D
chrY    57227301    57227400    Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D   Hsapiens-jaspar2016-TGIF2-MA0797.1  Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D
chrY    57227401    57227415    Mmusculus-HOCOMOCOv10-TFDP1_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-TEAD4_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-ZBTB6_MOUSE.H10MO.D

First, I converted all blank cells to 0. Now I'm wondering how I can replace the text in the cells with 1.
awk 'BEGIN { FS = OFS = "\t" } { for(i=1; i<=NF; i++) if($i ~ /^ *$/) $i = 0 }; 1' adrbvbeT_copy.tsv > adrbvbeT_copyN0.tsv

head adrbvbeT_copyN0.tsv
chr1    1   100 0   0   0
chr1    101 200 0   0   0
chr1    201 300 0   0   0
chr1    301 400 0   0   0
chr1    401 500 0   0   0
chr1    501 600 0   0   0
chr1    601 700 0   0   0
chr1    701 800 0   0   0
chr1    801 900 0   0   0
chr1    901 1000    0   0   0

tail adrbvbeT_copyN0.tsv
chrY    57226501    57226600    0   0   0
chrY    57226601    57226700    Mmusculus-jaspar2016-Tcf21-MA0832.1 0   0
chrY    57226701    57226800    Mmusculus-HOCOMOCOv10-TF7L2_MOUSE.H10MO.A   0   0
chrY    57226801    57226900    Mmusculus-HOCOMOCOv10-STAT1_MOUSE.H10MO.A   0   0
chrY    57226901    57227000    Mmusculus-HOCOMOCOv10-ZEP2_MOUSE.H10MO.D    0   0
chrY    57227001    57227100    NA-SwissRegulon-SMARCC2.SwissRegulon    Mmusculus-HOCOMOCOv10-PBX2_MOUSE.H10MO.C    0
chrY    57227101    57227200    Mmusculus-HOCOMOCOv10-ZIC3_MOUSE.H10MO.C    Mmusculus-HOCOMOCOv10-TBX5_MOUSE.H10MO.D    Mmusculus-HOCOMOCOv10-PIT1_MOUSE.H10MO.C
chrY    57227201    57227300    Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D    Mmusculus-UniPROBE-Tbp.UP00029  Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D
chrY    57227301    57227400    Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D   Hsapiens-jaspar2016-TGIF2-MA0797.1  Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D
chrY    57227401    57227415    Mmusculus-HOCOMOCOv10-TFDP1_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-TEAD4_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-ZBTB6_MOUSE.H10MO.D

You can use if (!$i) -- an empty string is false, and 0 is false. And if the field has the value zero, we're setting it to zero anyway.

Why do you need to do this? It seems like an XY problem

Not an expert, but :
awk 'BEGIN { FS = OFS = "\t" } { for(i=4; i<=6; i++) if($i ~ /^ *$/) {$i = 0} else {$i = 1} }; 1'

Since you have a very specific file format, this code starts from field 4, checking whether the fields after it are empty or not. I assume also that know that you file has 6 fields max exactly, so I "hardcoded" this knowledge.

Great it's working... Thanks a lot, @Sam.

source: Stack Exchange

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

sort: votes


text:
Converting a data frame to multiple lists
I have data frame each row is a pathway and it involved genes like

like
> head(h)
    Release.of.Cancer.Cell.Antigens.1  ARID1A    ATM  BIRC3   BLM BNIP3   BRCA1   BRCA2 BRD3  BRD4 BRIP1 CCNA1 CCND1
V1        Cancer Antigen Presentation    ARG2   ATF3    B2M BATF3   BLK   CCL22 CCL3/L1 CCL4  CCL5  CCR5  CD14  CD19
V2      T-cell Priming and Activation ADORA2A  BATF3   BTLA  CCL2 CCL22 CCL3/L1    CCL4 CCL5  CCR2  CCR5  CD1C   CD2
V3 Immune Cell Localization to Tumors     A2M ADGRE1 ANGPT2  APOE  ATF3     AXL   BATF3  BLK CCL13 CCL14 CCL18 CCL19

I want to convert that like
gene.sets <- list(Release.of.Cancer.Cell.Antigens = c("ARID1A", "ATM",  "BIRC3",    "BLM"),
            Cancer_Antigen_Presentation = c("ARG2","ATF3","B2M"))

I have not done anything so far as I am not certain what exactly I should google for.
By code the user kindly provided I get this
gene.sets = lapply(1:nrow(h), function(x) unname(unlist(h[x,-1])))
names(gene.sets) = h[,1]

[91] "IGF2R"     "IHH"       "IL11"      "IL11RA"    "IL12RB2"   "IL15"      "IL17A"     "IL1A"      "IL1B"
[100] "IL22RA1"   "IL2RA"     "IL2RB"     "IL6"       "IL7R"      "IRF1"      "IRF2"      "IRF3"      "IRF4"
[109] "IRF7"      "IRF9"      "ISG15"     "ITGAE"     "JAK1"      "JAK2"      "JAK3"      "KIT"       "LAG3"
[118] "LAIR1"     "LIF"       "LTB"       "MX1"       "NCAM1"     "NKG7"      "OAS1"      "OAS2"      "OAS3"
[127] "PDCD1"     "PDCD1LG2"  "PF4"       "PIK3CA"    "PIK3CG"    "PIK3R5"    "PRF1"      "PRLR"      "PROM1"
[136] "PSMB10"    "PSMB5"     "PSMB9"     "PTEN"      "REN"       "RORC"      "RSAD2"     "SELL"      "SELP"
[145] "SIRPA"     "SOX11"     "SOX2"      "SPRY4"     "STAT1"     "STAT2"     "STAT4"     "TBX21"     "TIGIT"
[154] "TNF"       "TNFAIP3"   "TNFRSF10B" "TNFRSF10C" "TNFRSF11A" "TNFRSF11B" "TNFRSF14"  "TNFRSF17"  "TNFRSF18"
[163] "TNFRSF1A"  "TNFRSF1B"  "TNFRSF4"   "TNFRSF8"   "TNFRSF9"   "TNFSF10"   "TNFSF12"   "TNFSF13"   "TNFSF13B"
[172] "TNFSF18"   "TNFSF4"    "TNFSF8"    "TPSAB1/B2" "TRAF1"     "TSLP"      "VEGFA"     "WNT7B"     ""
[208] ""          ""          ""          ""          ""          ""          ""          ""          ""
[217] ""          ""          ""          ""          ""          ""          ""          ""          ""
[226] ""          ""          ""          ""          ""          ""          ""          ""          ""
[235] ""          ""          ""          ""          ""          ""          ""          ""          ""
[244] ""          ""          ""          ""          ""          ""          ""          ""          ""
[253] ""          ""          ""          ""          ""          ""          ""          ""          ""

x=lapply(gene.sets, function(z){ z[!is.na(z) & z != ""]})

I don't know whether this will actually work since you didn't provide a reproducible example...
gene.sets = lapply(1:nrow(h), function(x) unname(unlist(h[x,-1])))
names(gene.sets) = h[,1]

Thanks worked. As different pathways have different number of genes involved to I got a lot of empty thing. The pathway with the most number of genes has 292 genes so the code assumes each pathway should have 292 genes that is why for the rest puts "". I have added the result of your code in my main post

Thanks again I have added solution to remove ""

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


text:
Loading scATAC-seq data containing matrix.mtx, fragments.tsv and barcodes.tsv files from GEO into R?
I want to analyse the single-cell ATAC-seq data from 10X downloaded from GSE158398. The files attached to this dataset are matrix.mtx.gz, fragments.tsv.gz, and barcodes.tsv.gz files, in addition to a peaks.bed.gz file.
How do I load this into R and make a Seurat object? Normally there is a matrix, barcodes and features files as output from scATAC-seq 10X data, so how do I handle only having a fragments file instead?
Thank you very much for your help with this!
Edit:
This is the code I would normally run to load the ATAC counts (matrix), barcodes and features and make the Seurat object. However, with this dataset, I do not have the features.tsv file, but a fragments.tsv file instead and so I am unsure how to change my code to accommodate this change in input data.
ATAC_dir <- "~/GSE169246_PDL1/ATAC/"

barcode_path_ATAC <- paste0(ATAC_dir, "barcodes.tsv.gz")
features_path_ATAC <- paste0(ATAC_dir, "features.tsv.gz")
matrix_path_ATAC <- paste0(ATAC_dir, "matrix.mtx.gz")

barcode_names_ATAC <- read.delim(barcode_path_ATAC,
                            header = FALSE,
                            stringsAsFactors = FALSE)$V1

feature_names_ATAC <- read.delim(features_path_ATAC,
                                 header = FALSE,
                                 stringsAsFactors = FALSE)$V1

mat_ATAC <- readMM(file = matrix_path_ATAC)

colnames(mat_ATAC) <- barcode_names_ATAC
rownames(mat_ATAC) <- feature_names_ATAC

chrom_assay <- CreateChromatinAssay(
  counts = mat_ATAC,
  sep = c(":", "-"),
  genome = 'hg38',
  fragments = "peaks.bed.gz",
  min.cells = 10,
  min.features = 200
)

seurat_ATAC <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks"
)

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

Welcome to the site Ebla. Could you kindly describe the code you have written so far to load the data?

Thank you, I have edited my original post now. I hope it is more clear now, please let me know if I can make it clearer


text:
Human genome similarities
As far as I know, we as humans share 99.9 similarity in our genomes. Further, we inherit 50% of the genome from our parents. So, it is correct to say that parent children have a 99.95% similarity in their genome?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Where is your source on the 99.9% number, and how do you know it's not 99.9134% or 99.8976% or 99.9999%? Also, what are you comparing - 99.9% shared with whom? A total genetic stranger?

The answer is no if the question is read at face value,

I don't know the metric used (is this just SNPs?), nor which part of the genome is included in the calculation, just the exome?
Moreover, if the genetic identity between parents is unknown, then the similarity of their offspring to either parent cannot be calculated.
Your calculation excludes de novo mutations, particular at neutral sites and copy number variation, e.g. microsatellites, however its not clear whether the data is simply the exome.
Homologous recombination is excluded, i.e. particularly occurring in the parental germ-line, there must be stochastic effects even in the first generation.
Finally, mtDNA is maternally inherited.

0.1% is the common rate of base substitution across the whole genome. 2) the number roughly varies between 0.07% and 1.2% across populations. 0.1% is considered a decent ballpark. 3) the de novo mutation rate is $10^5$ smaller and can be ignored, though CNVs etc may contribute a lot. 4) Recombinations don't matter. Exactly 50% of your autosome is inherited from one of your parent. 5) is correct, but sex chromosomes make bigger differences.

Thanks, you're saying the figure is singly SNPs at the chromosome level. "Exactly 50% of your autosome is inherited from one of your parent": we're talking homologous recombination in the parental germ-line, gamete formation. Homologous recombination is a major player in eukaryotic genetic diversity, unless the diversity of the markers used (SNPs) is so low it's difficult to track via linkage disequilibrium, i.e. breaks in physical linkage in a single generation. It will still be occurring. I dunno about homologous recombination in early mitotic events in embryogenesis.

Recombinations just switch the haplotypes to transmit. This doesn't change the calculation when you compare diploid genotypes.

It's a (over?) simplification. LD across a physical chromosome is a core calculation in eukaryotes. Humans too?

You are overthinking. I suggest you write down the calculation on a paper. It is simple. Recombination or LD is irrelevant to this question.

comment to (1),  according to https://www.nature.com/scitable/popular-students-page/132/ , coding region and non-coding region are included in the calculation of the 1% DNA sequence difference.

comment to (3),  in humans, a mutation rate of about 10 to the power -8 mutations/bp/generation was inferred from projects where both parents and their children were sequenced at high coverage.  (  http://book.bionumbers.org/what-is-the-mutation-rate-during-genome-replication/ )

The common wording is that the difference between two human haplotypes is 0.1%. This number only counts base substitutions, not insertions/deletions (indels) or structural variations (SVs). SVs together affects more bases than base substitutions. Also, the number is smaller is non-Africans (~0.07%) and higher in Africans (~0.12%). Anyway, let's denote this number by $\theta$, also known as the scaled mutation rate.
Because human is diploid, you need to count genotype differences when you compare two human individuals. Under the Wright-Fisher model, the rate of genotype difference on autosomes is
$$\theta\cdot\left(1+\frac{1}{2}\cdot\frac{1}{2}+\frac{1}{3}\right)=1.58\theta$$
The rate between a child and his/her parent is, I think:
$$\theta\cdot\left[1\cdot\frac{1}{2}+\frac{1}{2}\cdot\left(\frac{1}{3}+\frac{2}{3}\cdot\frac{1}{2}\right)+\frac{1}{3}\cdot\frac{1}{2}\right]=\theta$$
Note that I am sure the rate can be theoretically calculated but I am not confident in my calculation. Also, the rate of genotype difference on sex chromosomes or mitochondrial DNA is different.


text:
Check mutation of a gene

Hi all, I am checking if a gene which around 50k base pairs has mutation. So is that a variant from A to G? Would you please tell me how to do so? Thank you so much!
The reference sequence from GRCh38, sequenced using Illumina technology.

Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on [meta], or in [chat]. Comments continuing discussion may be removed.

Assuming this is a diploid genome, this looks like it could be a heterozygous variant (i.e. A is present on one chromosome, and G is present on the other chromosome).


text:
How do I calculate inbreeding on a single VCF file?
I would like to know the level of inbreeding for a single sample VCF file. I'm agnostic to the specific method of measuring inbreeding; any should work. However, I'd like to avoid having to specify a "background" or "reference" population. Ideally, the only input should be the single VCF file containing one human genome sample (WGS, hg38). (I will run this multiple times on different samples, but I'd like to run them individually.)
The existing tools seem to all require multiple VCFs.
Relatedly, is runs of homozygosity an acceptable metric for evaluating inbreeding in this context? And are there any special considerations to consider for various ethnic groups, like Ashkenazi Jews, that might influence the interpretation?

https://samtools.github.io/bcftools/howtos/roh-calling.html will output data, but I don't know how to get its output into a usable format. Maybe the autozygosity is just the number of lines in the ST section divided by the number of lines processed?

Possible solution:
for record in vcf_reader:
    if not record.is_snp or len(record.ALT) != 1:
        continue  # skip if not a biallelic SNP

    gt = record.genotype('Sample')['GT']

    if gt == '0/1' or gt == '1/0':
        obs_hets += 1
    elif gt == '0/0' or gt == '1/1':
        obs_hom += 1

num_variants = obs_hom + obs_hets
exp_hets = num_variants * 0.5  # under the assumption of Hardy-Weinberg Equilibrium

f = 1 - (obs_hets / exp_hets)  # calculate inbreeding coefficient

Example output: "-0.25"
Interpretation: We observed more heterozygosity than we expected.
I don't know if this method is valid, however.
A different way to calculate it:
for record in vcf_reader:
    if not record.is_snp or len(record.ALT) != 1:
        continue  # skip if not a biallelic SNP

    gt = record.genotype('Sample')['GT']

    if gt == '0/1' or gt == '1/0':
        obs_hets += 1
        allele_counts['A'] += 1
        allele_counts['a'] += 1
    elif gt == '0/0':
        allele_counts['A'] += 2
    elif gt == '1/1':
        allele_counts['a'] += 2

total_alleles = sum(allele_counts.values())

# calculate allele frequencies
p = allele_counts['A'] / total_alleles
q = allele_counts['a'] / total_alleles

# calculate expected heterozygosity
exp_hets = 2 * p * q * total_alleles / 2  # the total number of variants

f = 1 - (obs_hets / exp_hets)  # calculate inbreeding coefficient

Example output: "-0.45"
Interpretation: We observed more heterozygosity than we expected. I don't have a good idea of the meaning of this number.


text:
Missing samples when using function plotProfile() from DiffBind for ATAC-seq
I have 4 samples which are 4 conditions but when I made this plot, I got only 2 conditions.

Here is the code I used:
samples <- read.csv("path/to/sample_sheet_DiffBind.csv")
result <- dba(sampleSheet=samples)
result <- dba.blacklist(result)
result <- dba.count(result, bParallel=FALSE)
result <- dba.normalize(result)
result <- dba.contrast(result, minMembers=2)
result <- dba.analyze(result)
dba.plotProfile(profiles)

I followed DiffBind manual but not sure why. Thank you so much! This is the plot before analyzing.


text:
Good practices for doing sanity checks when parsing multiple PDB files
This is my first time posting to stack exchange, as well as to using data-driven approaches, so please correct me if my questions is not clear enough.
I am trying to parse 200+ .pdb files to extract residues within 15 angstroms of a certain ligand. I have managed to extract the resnames in a text file.
In order to do a sanity check, I first ran my code on the first 5 pdb files of this list and compared the outputs to expected output. However, I was wondering if there are any other steps I can take to make my check more robust.
Are there any 'gold standard'/good  practices related to this?

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

Welcome to the site. Could you kindly provide further details, e.g. are you using pymol, biopython PDB and what code have you written to date? Any reproducible example is always welcome.

Add a counter to be sure you're parsing the exact amount of pdbs you expect to parse ? Check for the ligand before extracting the interacting residue ? Check for duplicate pdbs ?

Thank you @M__. I am using VMD to parse through the PDB files. I can share my algorithm as follows:

Load the PDB into vmd. \

Identify the number of occurances of ligand, say Na+. \

For each of the occurance of Na+, identify the atoms within 5A of distance from it. \

Extract the resid and chain names for each of those atoms and store it as a pair in a variable pair_list. \

Filter out the repeated pairs to make sure I am counting every residue only once bold(quick Q: do you think the combo of resid and chain are sufficient to avoid repeated counting?) bold

Extract the resnames of those unique residues and save it in a text file.

Repeat this process for all the PDB in current directory.

I have added the TCL script doing this for NH4+ at this github link: https://github.com/priyanshurg/TCLScripts/blob/main/reslist4-Within5.tcl

Thank you @pippo1980ONSTRIKE!
Could you recommend what property in the PDB file I can use to classify PDBs as similar? For example there are many versions of the KcsA protein in RCSB.
Your point about duplicate pdbs made me realise that my database may be biased by proteins that are studied more often compared to others. Do you have any thoughts on how you'd go about mitigating this bias?

dont know see  https://actanaturae.ru/2075-8251/article/view/10699/pdf materials % methods, you could drop structures too similar or sequences too similar I believe its on you , try both and look at the differences  or you could just concentrate on the active site alone using similar approach


text:
How to determine the phylogeny of a specific set of genes
My goal is to understand the phylogeny of genes involved in bacterial restriction-modification systems, in psychrophilic bacteria.
To do so here is what I have done.

Downloaded the protein sequences of all RM genes in psychrophilic bacteria from REBASE. I also downloaded a set of genomes I'd like understand the phylogeny of from NCBI.
Clustered these genes using usearch's uclust method.
Aligned the resulting centroids using MUSCLE5.
Built an HMM from the aligned centroids using HMMER3.
Ran hmmsearch to get hits from my HMM on each genome (AA).

I am unsure how to proceed now. I now a phylogenetic tree requires an alignment, but I'm not sure where my HMM hits fit-in. Am I meant to use hmmalign on the output of hmmsearch and then feed the resulting alignment to IQ-TREE or similar?
Thank you for your help as I learn bioinformatics.
PS: Reading through the forum it seems like MOLE-BLAST does exactly this, however I have too many genes and their length is too long to use the MOLE-BLAST as hosted by NCBI at least.
I should mention I set the identity threshold of search to 50% for proteins. Also, I got this method from this paper which achieved something similar (https://pubmed.ncbi.nlm.nih.gov/30893937/).

Your strategy is fine but the order is wrong.
What I would do is:

2 -> checked by 4 and 6

Run usearch and check the output of that in HMMER3. You can then use this to kick out stuff from usearch, that are not real homologues.  Thats what HMMER3, hmmeralign (I dunno hmmersearch), are doing (below). It's also possible the HMMER approaches could replace usearch.
Once you have identified your homologues with confidence you can then align each separate group via muscle5
Your way
What you are doing is using HMMER3 to assess the quality of the alignment. What could be a problem is feeding usearch straight into muscle5 ... a non-homologous protein will get you a bad alignment. The bad alignment will not give you a good tree. Both aligners and trees need nice homologues, otherwise they don't work. The errors between aligner (muscle5) and tree are also synergistic. Thus if that bad sequences isn't removed by hmmer stuff then alignment errors compound onto sequence errors in the bad sequence giving a really bad tree.
You could of course align, HMMER3, filter stuff then align again (thats much better a bit complicated).
I don't know how good usearch is - but you are relying on it alot, supposing it gives you orthologues at <25% identity that could be any protein. Basically you are feeding that stuff straight into muscle5 and if HMMER3 doesn't trap it that will go straight onto a tree. That could be a bad alignment and that would be a bad tree.
Summary

What you have to understand is phylogenetic trees are very sensitive to bad sequences, i.e. stuff that is not true homologues. The aligner will also be sensitive. There will also be bad synergistic relationships between the two.

HMMER3 and is associated methods (hmmeralign) are trying to clean stuff up before that stage, rather than simply relying on usearch. Prof Eddie (HMMER and pfam[?] author) has an excellent reputation in this area - not in tree building but identifying protein families and thats what you are doing here.

hmmeralign doesn't replace muscle5, but it will provide a very good assessment of the protein families around a given seed. It's going to give you the homologues. That output would then go to muscle5 to get homologous sites.

The tree needs both homologues and homologous sites.

The whole hmmer project has a really good reputation in the area of protein family identification and thats why its been recommended.
Probable simple answer I think the answer to the question is you're are confusing hmmeralign with muscle5. They are complimentary but distinct. Again hmmeralign identifies homologues and the muscle5 designates homologous sites.

Answering the comments

50% cut-off

It looks like usearch clusters according to a %age identity threshold.
There are three weaknesses:

Gene duplication can take place about 50% ... it depends whether thats cool, usually it isn't, but it's your call.
There can be complete change of function >50% identity and this is about the definition of a homologue. If you are using a functional definition its not a homologue
Homologues can have <<50% identity and thats very common certainly in viruses I don't know about bacteria.

hmmersearch is going to rank your data according to the hmmer algorithms. It's going to be like blast,but all hmmer algorithms are focused around defining protein families (it's a better approach than blast for this application). So it provides another option to assess the output of usearch. Thus if there a gene duplication for example you should be able to spot it (it will be a sharp break in the continuous output). This will help you assess you 50% cut-off.
hmmeralign would be quite cool, because you can take a representation of usearch and then search that via hmmeralign and see if you get a similar output. It is extremely important for assessing 50% cut-off, does the homology drop below 50%? hmmeralign will tell you, but it must see the complete original database (not the stuff extracted from usearch. You can cross check that with hmmersearch to assess where hmmeralign is making the cut-off. So hmmeralign will independently (or usearch) will find you homologues and hmmersearch will give you an insight about how it's done it.

I have trouble seeing where hmmsearch fits in here, or is it not meant to?

One more point of confusion. If you suggest I check hmmalign with usearch, then what alignment do you recommend I feed to hmmbuild? Just a straight alignment from muscle5?

I would go straight into muscle5, without hesitation. Thats because if something goes wrong I'd spot it in the tree and the outputs. If you're phylogenetics isn't so strong, you might want to perform all 'sanity checks' in hmmer.


text:
Mean and median identical with no variation in FASTQC results
Mean and median identical with no variation in FASTQC results
After using FASTQC, in the section of "Per Base Sequence Quality" I encountered a problem when I saw the Box & Whisker type plot. The median value and the mean quality are the same and I have no yellow box, as shown here:

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

Can you show the first few records from your fastq file and say roughly how many reads you have? Where did the fastq come from? Is it possible that 'dummy' quality scores of Q30 were added by someone?

Please look to comply with @Chris_Rands. I don't think anyone will answer until this request is fulfilled.


text:
Pal2nal translation of large multi-fasta files produces a codon translated file where some sequences half length of the average
I did sequence alignment of a large peptide multi-fasta (n= 4991 sequences). The peptide alignment has sequences with the same length and pal2nal went through just fine... except some of the codon sequences are at half length. If average is X then some sequences are X/2. This is choking IQ-Tree.
I have tried both MUSCLE super5 and MAFFT. The error remain the same (i.e. MUSCLE or MAFFT both lead to some sequences having half of average length) except for different average lengths and average half length in MUSCLE and MAFFT codon sequences. I have pulled out and played with the sequences causing the issue and they seem to be in frame.
Example of peptide sequence not causing an issue: RKVEAFLLFKEMGERGCQPNVHTYTVLIDSFCKERNLDDARKLFDDMFKKGLVPSVVTYNALIDGYCKEGMTEAALEILGMMESKKCNPNARTYNELICGFCKAK
corresponding cds
AGGAAAGTGGAAGCTTTTCTACTTTTTAAAGAAATGGGTGAAAGAGGTTGTCAGCCTAATGTTCATACATACACTGTGCTTATTGATTCCTTCTGTAAGGAAAGGAATCTTGATGATGCCAGGAAATTGTTTGATGACATGTTTAAGAAAGGTTTGGTTCCCAGTGTGGTCACTTATAATGCTTTAATTGATGGGTATTGTAAAGAGGGAATGACTGAAGCTGCATTAGAAATTTTAGGTATGATGGAATCAAAGAAATGCAACCCTAATGCTCGGACCTACAATGAATTGATCTGTGGATTTTGTAAAGCTAAA'
Example of peptide causing issue:
GLCKGGRLNDAWEIFQYLLAKGYQLNVHTYNAMVHGFCKEGLLDEAISLLYKMEENGCVPNSVTFNVVL
corresponding cds
GGTTTGTGCAAAGGTGGTAGATTAAATGATGCGTGGGAGATTTTTCAGTATCTTTTAGCGAAAGGTTATCAACTAAATGTCCATACATATAATGCGATGGTTCATGGTTTTTGCAAAGAAGGTTTGCTTGATGAAGCAATCTCCCTGCTTTATAAAATGGAAGAGAATGGTTGTGTCCCTAATTCTGTAACTTTTAATGTAGTCCTT
Any idea what might be going on?

The example looks fine to me. The peptide causing an issue is 69 residues in length and the cds is 207, which is correct. If the peptide is supposed to be ~140 residues we'd need to see the full length peptide and/or nucleotide. BTW what is the proportion of half lengths against full transcripts?

The peptides are equal in length once they are aligned alongside their orthologs/homologs (i.e. on the other side of sequence alignment). With MAFFT the peptides become something like 31k (including gaps) characters long. Post translation the sequences are something like 93k long and the ones out of whack are about 54k (all of them, the same exact length). Sorry for the confusion, the peptides are not supposed to be the same length per-alignment.

I can see whats happening that the 'issue' sequence is a locus within the full gene for Acacia spp.. The thing is I can't see the original sequences, i.e. before Pal2Nal. The reality I can probably solve this with all the information. At present I can make a guess, but it's no better than a guess.

I am happy to provide the original data, it's all open-source. Do you have a preferred way to receive them?

I don't want to wade through 5K sequences. What I am really asking for is a reproducible example. The example provided doesn't help understand the issue, it's the input into pal2nal I'm looking for. What I want to do is compare the input vs the output.

The first 8 sequences, containing both 'working' vs 'issue' sequences are here. It's a fork of Firefox send, I can use something else if you are uncomfortable with that. By locus if you mean alternate splicing then it is possible that that is what is at play.

Not a problem, doing global alignment on a phylogeny that spans a billion-years was a gamble to begin with. I am grateful you took the time to look at it. I suspect these homo-logs are not real or if they are then the global alignment algorithms we have might not be ready to work on something a phylogeny this deep.

I had a quick blast of the example sequences and this is the Acacia so ancient vascular plant. It forms a locus within the protein in the first part of the question. This protein is comprises a large repeat amino acid motif (check NCBI), the second "issue sequence" described in the question was close to the size of the repeat locus, but start/stop of the locus is not in accordance with the repeat structure. Thus I don't think the repeat locus is part of the issue.
Genome assemblies I noted the examples you gave were all draft genomes and this suggest the problem was with the original genome assemblies. If this was true the way to overcome this would be to reassemble the sequence using a reference genome, bwa-mem2 type thing. Thus, here the error is occurring way before the sequences is even downloaded. Just to note, repeat structures can often be collapsed on a draft genome assembly, but I don't think thats happened here.
Reading frame what I think you are doing is performing an alignment in the proteins and then using that alignment as the template which is then transferred to the sequence data. So now the nucleotide data is nicely aligned including all indels. That appears to be the goal of pal2nal. A common error here is reading frame error, particularly at the beginning of the sequence. I use transalign regularly for this purpose and do encounter that error. I don't know how pal2nal deals with reading frame errors so I can't comment, but thats why I wanted to see the nucleotide input to pal2nal to check this.
Trees the reason this is taking place is because iqtree doesn't permit a tree. That I can help with, I know about tree building. If you post your iqtree command as a separate question , together with the iqtree error, I can probably answer that. What I'll identify is how to build a tree, this likely to exclude sequence, or could be addressed as two trees. Dunno, but if the bottleneck is can't build a tree - thats easily fixable. In this case, is the sequence issue an assembly error, reading frame error ... something else. There's not enough data that I could point to as the underlying source.


text:
How do you extract the full list of files in a GEO archive?
I am new to use of GEO and having trouble figuring out a straightforward way to grab the names of the files,  (through programmable handles like GEOparse or otherwise) that is the entire and complete list of all files fastq and otherwise (NOT downloading the actual data, but the names and descriptions) within a specific GEO archive (GSEXXX). Albeit with the intention of categorizing the data into  more understandable and friendly to browsing segments, the data presentation in a GEO archive seems very fragmented, making it hard to find what you are looking for. When you submit your data to GEO, they require you to have a spreadsheet, that is essentially the list of files you are submitting with their corresponding meta-data descriptions. Is it possible to get that exact submission spreadsheet associated with a data archive?

Thank you very much @Sam, that looks very relevant. Just a superficial exam and experiment with geofetch has proven more helpful than GEOParse approach. I am not able to get the list of fastq files using --just-metadata so far, BUT I am able to  the processed file list such as peak files. I will play around with it more to see if I can crack fastq files also and I think that would completely answer my question.

Unlikely that you will be able to get the original submission spreadsheet, but you can use Geofetch


text:
Different human reference genomes from NCBI and ENSEMBL, and variant annotation with Variant Effect Predictor (VEP)
Background: I have been working on a human variant calling pipeline, from whole exome reads to variant annotation.
For variant annotation I used the Variant Effect predictor (VEP), which is Ensembl based.
I have used the NCBI reference genome GRCh38.p14 and the older version GRCh37.p13. The latter version was for comparison with older variant calling pipelines. It is about the latter geneome version I want to continue the discussion.
Problem: I sparked a little discussion due to the use of different versions of "hg19" genome (GRCH19, B37, NCBI, Ensembl) in a research project. This is related to this article. Due to discrepancies in these versions, I also wholeheartedly believe these cannot be used interchangeably. But this also got me thinking about using VEP, based on Ensembl and likely their reference genome, while I also used in my pipeline the reference genome of NCBI. When I look at just two example genes, you can already notice the difference:

Question 1: is VEP able to correctly annotate variants likely based on the Ensembl reference genome, while I am using the NCBI reference genome?
Question 2: any thought/knowledge on using the one or the other reference genome, Ensembl or NCBI?
Question 3: Follow-up on Question 2, how do these versions relate to the UCSC version? There does not seem to be a cache for this version?

VEP asks you to pick a specific cache when you annotate, so you choose if you want NCBI or EnsEMBL. I prefer the latter because it is comprehensive

Hi, please don't answer the question in the question body, write it as a seperate answer to the question.


text:
Louvain/Leiden community clustering implementation in C/C++ for gene expression data
Does anybody know of a C/C++ based implementation of the Louvain or Leiden algorithms that will take gene expression looking csv input files, with the genes in columns and the samples in rows, and perform community clustering on it? Something like
@lh3lh3 https://github.com/lh3/bhtsne but for the community clustering algorithms.

igraph does Leiden and Louvain clustering. It is a C/C++ library and used in R or Python.

Thanks, I looked at the igraph library but couldn't find out how to turn the matrix of input data (e.g. the 28x28 pixel data from mnist) into an input graph to run the community clustering method.

I see. Personally I wouldn't use unsupervised learning for images (e.g. the MNIST data set), I'd use supervised learning: XGBoost is really cool for this BTW. Thus I dunno.


text:
Keeping DNA sequence after changing FASTA header on command line
I have a FASTA header that looks like this:
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2 runid=606f271fc97598006ba5a922136a2c304cef75a5 sampleid=Pool12-1 read=19008 ch=301 start_time=2021-07-03T08:48:18Z barcode=barcode01

And I am able to change it to the desired output here:
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2_001

Using this command:
grep '^>' 001_old.fasta | cut -d ' ' -f 1,8 | sed '/^>/s/$/_001/' > 001_new.fasta

However it completely deletes the DNA sequence for that read and only outputs the names of the headers to the new fasta file:
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2_001
>815440d5-e033-446c-b7dd-5fad4105184a_001
>d36a0fea-101d-4bd9-bee4-adec7fe4d43c_001
>38ca9eb8-bbd2-4d0a-ab06-79d134ad46f1_001
...

I think it has something to do with the cut command I'm using. I tried just doing 1 but that gave the same output so I figured if I included the 8th position of text after a space it would work but again no dice. If someone could help me to include the DNA sequence under the FASTA header in this new file that would be great.

Do yourself a favor and use Bipython for parsing and modifying FASTA files.

The grep, cut, sed pipeline could also be replaced with the following awk command:
awk '/^>/ { print $1 "_001"; next }1' file.fa

The general syntax of an AWK command is condition { action } where condition specifies a condition or regular expression, and action specifies a block of code or set of commands to be run when the condition is met. Both condition and action are optional. Here, we use a regular expression to match lines beginning with a ">" symbol, and print the first whitespace-separated field followed by "_001". We then use next to skip the remaining actions and move to the next line, which avoids any further processing of the current line. The 1 on the end is actually another condition that always evaluates to true. Since there's no action associated with this condition, AWK performs the default action which is to print the entire line. This effectively prints all lines that do not begin with a ">" symbol.

thank you this worked perfectly! and thank you for the explanation too

The problem is in the grep command you're using, which only sends the fasta file headers to the output. Here's my guess at what the commands in your pipeline do:

grep '^>' 001_old.fasta - extract lines starting with '>' from 001_old.fasta

cut -d ' ' -f 1,8 - assume file contains space-delimited records; extract only the first and eighth field from each line

sed '/^>/s/$/_001/' - in lines that begin with a '>', replace the end of the line with '_001'

> 001_new.fasta - write result out to 001_new.fasta

Of all of these, I think only the sed and output parts are important, but it needs to be slightly modified to work as intended:
sed '/^>/s/ .*$/_001/' 001_old.fasta > 001_new.fasta

Assuming I've got the syntax right, this should search for lines beginning with '>', and replace everything from the first space up to the end of the line with '_001'.
I'm more familiar with perl, so here's a perl one-liner that should do the same thing:
perl -pe '/^>/ && s/ .*$/_001/' 001_old.fasta > 001_new.fasta

thank you so much! this makes a lot of sense to me as well.


text:
Correlation heatmap of RNA-seq clusters all samples together leading to very low no. of DEGs
I am writing to you to take an input or may be you can provide a different perspective. I am at wits end :(
So I have nearly 200 samples. These are separated into two groups (group I treated with drug1 & group II treated with drug II) of course these two groups have sub groups of responders to drug and non responders.
I have used DESeq2 to analyze this data. The first thing that I notice about this data is that the in the correlation plot the samples clustering across the groups. In past I have seen some samples clustering across the groups but this a lot. There is no cluster specific to one group.

Secondly, I get like 5 or 6 DEG that are statistically significant and have fold change of 2.
I have been part of similar studies with same tissue and comparable drugs but have never noticed something like this.
I am not sure what is happening here and why it is happening.
Any help or insights are sincerely appreciated.

What type of samples are these - Are these human samples? Cell lines? Different tissues? Different animals? Have you done any filtering or correction for whole-genome genetic similarity (i.e. ancestral similarities)?

This answer has two stages,

Diagnostics, what analytically has gone wrong?
Where to go from here.

I am pretty confident about point 1. I've been in comparable situations and I know trees. Stage 2, is more open to opinion, however there's no point criticising something (point 1) without pointing out how it can be corrected (point 2).
If your question is within group heterogeneity, excepting for the large group of non-responders associated with D1, thats a different answer.
Low number of DEGs is about your control group? Difficult to comment because the DEG calculation is not declared, nor volcano plots, cut-offs, FDR etc...

You have to keep in mind this data mining and the there are some interesting groups emerging. For first stage analytics looks reasonable (low DEGs dunno).
Result You have 9-10 clusters and those are informative.  You can see this because they are forming a checker-board across the heat-map.
Issue The issue is that the deeper branches below the 9-10 clusters are not informative. Do see? The deep branches are likely to be repositioned in different orders in different analyses. Those alternative order might make better biological sense. On the heat-map shown there's a risk of attaching two much important to deep branching orders for a UPMGA tree for a correlation dendrogram.

Heterogeneity across the matrix is not permitted in UPGMA and causes artefacts - because you are trying to say there's a uniform rate when in most biology that assumption doesn't work. This will usually result in incorrect branching orders.
You are forcing the data into a bifurcating tree: it might not be at all, could be a large polytomy (star tree).
Even data sets designed for trees (correlation matrixes are not good data sets for trees) ... branching orders go wrong for loads of reasons.

The second stage analysis would be to take the correlation matrix and apply the neighbour-joining (nj) algorithm to it and see whether the deeper branch align differently. There are also much better tree methods, but that starts getting complicated fast, because there are alternative approaches, such as co-variance. Advanced tree methods (Bayes and ML) will model shifts heterogeneity - however it needs some expertise to get it working.
I'd start assessing the assumptions UPGMA vs neighbour-joining does that take me closer to a biological explanation, then a 9x9 checkboard? Thus would the number of squares in a checkerboard by greatly reduced via nj?
Machine learning Supervised learning (ML), is the trendy way into the data set and in this day an age thats probably the way I'd go. What I'd do is use the 9-10 clusters here as an initial assessment of the training targets. I'd then perform whats called a "feature selection analysis" and assess whether that output makes biological sense. The feature selection should start identifying what parts of your data are causing the shifts between those 9-10 clusters.


text:
How do I add reference/alternate alleles to my vcf file?
Here are the first few columns of my VCF v4.2 file. As you can see, the alt allele is missing. The VCF's refs and alts don't correspond to the reference genome.

#CHROM
POS
ID
REF
ALT

9
77051
rs9408670
T
C

9
114350
rs12685480
A
.

9
116581
rs558835286
G
.

9
117550
rs267602077
A
.

9
117720
rs10959293
T
.

9
126903
rs62533602
G
A

9
133828
rs78804900
C
T

9
154736
rs200683662
T
.

9
154933
rs17795008
T
G

9
174425
rs2998944
T
C

9
174500
rs142581442
G
.

9
175050
rs62532020
G
A

9
185632
rs2992854
C
A

9
188894
rs3951559
G
A

9
191843
rs3008128
G
A

9
202061
rs117675283
T
.

9
204201
rs10964134
G
A

The missing alt allele happens when the genotype is 0/0. The VCF file isn't saying anything incorrect. However, there is a tool I'm trying to use that filters out "0/."genotypes.
A better way of formatting them (which I'm trying to achieve) would have the ref/alt alleles correspond to the reference genome, and just have a redundant letter in the case of a homozygous genotype.
I'd like my VCF to display the correct alt/ref for that position even if there are no samples with a genotype corresponding to the alt allele. How do I do this? I'm not sure what this process is called.
Ideal output:

#CHROM
POS
ID
REF
ALT

9
77051
rs9408670
T
C

9
114350
rs12685480
A
G

9
116581
rs558835286
G
A

Are you sure you need to update the alt? The problem with 0/. genotypes is that period usually means “missing” so that genotype is kinda bizarre. There’s one allele we know, the reference, but also one we couldn’t figure out? If you think these should actually be homozygous reference genotypes, convert them to “0/0”.

If you mean the tool can’t handle “.” in the ALT column, maybe just use awk to update all the dots to “N”? That’s a valid VCF allele but won’t be a ref.

It would still filter that out, unfortunately. My solution for now is to create valid letters (ATGC) for the ref, replacing dots.

You should include the genotype column of the VCF, so that we have a clearer picture. As has been mentioned, it is not the same a missing (./.) genotype than a homozigous reference. The former representation usually means the actual genotype could not be establish by the variant caller and assuming hom ref might cause mistakes in dowstream analyses a

Let me know if you see any potential issues with this solution. It just replaces the dots with random letters (easier/faster than aligning to reference genome).
def add_ref():
    # Directory path
    dir_path = '/path/'

    # Function to process genotypes in each row
    def process_genotypes(row):
        genotypes = np.array([x.split(":")[0] for x in row[9:]])
        if row[4] == "." and not any(val in ["1/1", "0/1", "1/0"] for val in genotypes): # to make sure I'm not introducing an allele that isn't present
            # also can do if row[3] in ["A", "T"]: then random.choice(["C", "G"])
            if row[3] == "A":
                row[4] = "C"
            elif row[3] == "T":
                row[4] = "G"
            elif row[3] == "C":
                row[4] = random.choice(["T"]) # I know
            elif row[3] == "G":
                row[4] = random.choice(["A"])
        return row

    # Get all .vcf.gz files in the directory
    files = glob.glob(dir_path + '*.vcf.gz')

    for file in files:
        print(f"Processing {file}...")
        # Name of the new file
        new_file = os.path.join(dir_path, 'new_' + os.path.basename(file))
        with gzip.open(file, 'rt') as in_f, gzip.open(new_file, 'wt') as out_f:
            for line in in_f:
                if line.startswith('#'):
                    out_f.write(line)
                else:
                    row = line.strip().split('\t')
                    new_row = process_genotypes(row)
                    out_f.write('\t'.join(new_row) + '\n')

Okay I see. Python pandas is an alternative. If this approach runs out of memory its worth thinking about.

You best make sure to that the downstream program doesn’t use the alt allele field for something, or else it could cause problems

One of the downstream programs is imputation. Not sure if it uses the alt, but I don't see why it would.


text:
How can I give different names to files in a directory with a for loop in a bash script?
I'm expecting to get 17 different paired-end fastq files (34 in total), so I want to make a bash script to just run my code through all the fastq files in a directory at once. How can I change the name of the input and output files each time the script runs through each file? So when it moves to the file_002, all names have file_002 at the beginning instead of file_001, and so on. And also, when merging the R1 and R2 reads how can I make that it only merges the correspondant files with a loop? for examples merging only file_001_R1 with file_001_R2, file_002_R1 with file_002_R2, file_003_R1 with file_003_R2, and so on.
This is my code:
for file in directory_name
do
pear -f file_001_R1.fastq.gz -r file_001_R2.fastq.gz -o file_001.fastq
cutadapt -g TGATAACAATTGGAGCAGCCTC...GGATCGACCAAGAACCAGCA -o file_001_barcode.fastq file_001.fastq
cutadapt -g GTGTACAAATAATTGTCAAC...CTGTCTCTTATACACATCTC -o file_001_UMI.fastq file_001.fastq
seqkit concat file_001_barcode.fastq file_001_UMI.fastq > file_001_concatenation.fastq
seqkit rmdup -s file_001_concatenation.fastq -o file_001_unique_pairs.fastq
seqkit subseq -r file_001_unique_pairs.fastq > file_001_unique_barcodes.fasta
bowtie -q --suppress 1,2,4,6,7,8 -x ref_index file_001_unique_barcodes.fasta > file_001_barcodes_allignment.bowtie
sort file_001_barcodes_allignment.bowtie | uniq -c > file_001_barcode_counts.txt
awk 'BEGIN{print "Barcode,TF_variant,Code"}{print $3","$2","$1}' file_001_barcode_counts.txt > file_001_barcode_counts.csv
done

Thank you.

I assume the code is directory_name/*, $file is now by each file within directory_name. At the moment $file is simply directory_name as its written. There'll be bash experts here that will describe the loop for the paired end files.

Hi, please consider upvoting and accepting one of the answers. These are good answers from two different perspectives: one pulling the directory file contents; the second using a numerical loop. Both would work.

Here's what I use - along with an example you can modify as needed:
FASTQ_DIR="/directory/containing/FASTQs/"

for R1_FASTQ in "$FASTQ_DIR"*R1*.fastq.gz
do
                         BASE=$(basename "$R1_FASTQ" _R1.fastq.gz)
                         R2_FASTQ="${R1_FASTQ/R1/R2}"
                         R3_FASTQ="${R1_FASTQ/R1/R3}"
                         echo "BASE: $BASE"
                         echo "R1: $R1_FASTQ"
                         echo "R2: $R2_FASTQ"
                         echo "R3: $R3_FASTQ"
                         # process R1/R2 here using "$R1_FASTQ" "$R2_FASTQ" for input, and "${BASE}_XXXXXXX.fastq" for output
done

Basically, the strategy is to iterate over all the R1 FASTQs in a directory, calculate a BASE name that can be used to rename output files with arbitrary extensions (like .sam, .bam, etc...), and use BASH shell parameter expansion to calculate and set variables to the R2 (and R3 if necessary) reads present in an Illumina run that uses standard Illumina naming conventions.
This may need modification depending on the naming conventions you are using.

Something like this? This uses printf to create a zero-padded 3-digit number corresponding to 17 numbers from 1 to 17:
for fileBase in $(printf 'file_%03d\n' $(seq 1 17));
   do echo ${fileBase}_R1.fq.gz ${fileBase}_R2.fq.gz;
   pear -f ${fileBase}_R1.fastq.gz -r ${fileBase}_R2.fastq.gz -o ${fileBase}.fastq
   ...
done


text:
LD clump GRCh38 GWAS results
The vignette of R package ieugwasr describes a plink based wrapper function for LD clumping GWAS data using the 1000 genomes reference dataset (http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz). I would like LD clump GRCh38 based GWAS data. I am under the impression that the previous LD reference dataset is based on GRCh37. Therefore, I kindly ask:

Is 1kg.v3.tgz based on GRCh37?
Can you use GRCh37 based reference to LD clump GRCh38 GWAS results?
Is there a GRCh38 reference version of 1000 genomes that is split to different super populations (EUR etc)?

Is 1kg.v3.tgz based on GRCh37?

Is there a GRCh38 reference version of 1000 genomes that is split to different super populations (EUR etc)?

Yes, the genome build for 1000 Genomes phase 3 release is in GRCh37. Luckily ISGR has reanalyzed their original sequencing data based on GRCh38. The data is provided both at the ISGR website https://www.internationalgenome.org/data-portal/data-collection/30x-grch38 and it is also made available by Christopher at PLINK2 website https://www.cog-genomics.org/plink/2.0/resources. However I suspect the data was not split into super populations.

Can you use GRCh37 based reference to LD clump GRCh38 GWAS results?

I wouldn't recommend doing so. Some variants identified in GRCh37 has been merged into others in GRCh38. Some variants are assigned different rsIDs. The position for the variant may also differ in different genome builds. Therefore clumping GRCh38 GWAS results with GRCh37 based reference may lead to fewer variants correctly mapped, which may influence your furthur analysis.
Hope you find it helpful. Best regards.


text:
Gene Expression Database for Psychological Disorders
Is there a comprehensive database available online of gene expression across the brain taken from individuals with psychological conditions? Resources like the Allen Institute's Human Brain Map (link here for reference) feature heatmaps of gene expression taken from healthy volunteers, and I'm looking for similar datasets from patients with particular conditions, especially sensory processing disorders such as ASD and ADHD. Any and all data is appreciated!

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Welcome to the site. To be honest I would suggest you are likely best placed to answer this question, given your expertise demonstrated on other SE sites. NCBI SRA will host the datasets required and their reference number will be given in the associated publications. This is a standard approach to many data sets associated with RNAseq data, i.e. they are exclusively on the established generic repos.


text:
Handling duplicates in untargeted lipidomics
I have a MS-lipidomics dataset - from human samples split into a treatment group (n = 2) and control (n = 2) - which contains several detected targets that matched to the same lipid. I don't have access to any of the raw data, only the Tier 1 and Tier 2-matched lipid names and ID numbers.
There are around 300 duplicates out of a total of 1000 detected markers, so I'm reluctant to exclude them entirely, but I couldn't find any descriptions of best-practices for what to do in this situation. I'm considering testing the correlation between duplicates and removing the ones lie below some threshold (e.g: 0.5), then averaging the remaining duplicates - but I'm not sure if this would be a valid way of approaching this.
Are there any good ways of addressing these duplicate values?

Sorry how is the threshold derived thats quite key here?

I don't know lipids, but if 1/3 of your data set comprises duplication thats a skewed distribution with a low mean and high variance (of the duplications). It is likely to follow a negative binomial distribution. You need to plot a histogram of your data and you will find that 0 (no duplicates) is massively represented, 1 duplicate is next in frequency, 2 duplicates are lower frequency ... and on until a very small proportion of the data comprises a large number of duplicates. This is a very common distribution throughout biology.
The advice is:

Do not remove the duplicates because they are a vital part of the distribution, i.e. biological information.
Construct a histogram of each data set.
Don't take "an average", i.e. don't take the mean of skewed data, thats very important if it has a high variance.

Easy solution The easiest way to deal with this data is take a median and construct upper and lower quartiles around the median (IQR), especially the frequency (i.e. the distribution duplicates). This will give an approximate comparison between data sets. This alone might work and certainly gives you "cut-offs".
If it is a negative binomial distribution then testing the significance between two negative binomial distributions beyond the IQR is an advanced area of statistics and a separate question.


text:
What would be the best method to obtain every prokaryotic psychrophile genome?
I am interested in running an analysis on the genome of every psychrophilic bacteria. Here defined as bacteria that have an optimal growth temperature at less than 15°C.
What would be the best way to go about collecting these?
I have yet to find a "growth temperature" filter on databases like NCBI. Though REBASE does include that (however obtaining the genomes from that site does not seem straightforward).
Thank you

I ended up copying the GenBank (primary) column from the REBASE search results into an excel sheet. From there, I created a file with one accession per line called Psych_Accessions.txt. I then ran the following command using NCBI's command line tools.
epost -db nuccore -input Psych_Accessions.txt -format acc | efetch -format fasta > psych_genomes.fasta

The database to use is BactoTraits. This will help either in the entirety or at least provide a solid starting point. Its here
BactoTraits – A functional trait database to evaluate how natural and man-induced changes influence the assembly of bacterial communities
by Aurélie Cébron (2021) Ecological Indicators 130

We have created a user-friendly database of bacterial functional traits, thanks to the properties of 19,455 bacterial strains. This database has been called BactoTraits. For example, oxygen preference, size and shape of bacteria, motility, optimum and range of pH and temperature, genome GC percent and trophic type are among the 19 traits included in BactoTraits. Based on the best-informed strains in the database, we identified five functional groups (i.e. mesophiles, competitors, colonizers, stress-sensitives and stress-tolerants) exhibiting a wide strain taxonomic diversity but with quite similar trait profile combinations.

Thus it's nearly 20000 species so your batteria should be there.
You can download the exact database here. This is Univ de Lorraine's site, France. Its 9.4MB and called

BactoTraits_databaseV2_Jun2022.csv

The supplementary information in the paper gives all the details for interpreting the database columns.
Given your skills in development you should have no trouble converting the csv format into a functional database and interrogating/querying it. My understanding of the paper is that the authors are microbiologists rather than coders, which in context is perfect.

Linking NCBI
BactoTraits contains two strain designations

DSM no. (https://www.dsmz.de/collection/catalogue/microorganisms/catalogue)
ATCC no. for example  ATCC BAA 1260 is here on NCBI Granulibacter bethesdensis, which checks out on BactoTraits and is DSM 17861
txid364410[Organism:exp]  is the SRA number from the above link ... so its all there, and here https://www.ncbi.nlm.nih.gov/sra/?term=txid364410[Organism:exp]

So ATCC links to NCBI links to the SRA and genome.
DSM are not necessarily available on NCBI, but ATCC will there. NCBI do have lots of databases for download, a database comprising ATCC to NCBI genome ... dunno, but Entrez will be able to link the ATCC (American Type Culture Collection)

While this is a great database for getting the names of psychrophiles, I was hoping for something that would tie temperature optima to NCBI project numbers or other direct genomic artifacts to make it easier to download. I will keep this question open for a bit longer before marking your answer.

Comment answered about below the line @donkey


text:
building custom gtf file with entrez ID for any organism
How do I generate a custom gtf file format for  a organism of my interest similar to human gtf file, in the custom gtf file I want to add entrez id to the new custom gtf.
My organism is canis_familaris
I  have to use is this source.
How do I start, is there any tools there to build custom gtf file from the GENE_INFO file?
From the gene_info file I extracted the ensembl and the entrezID column so here as an example I have taken only the top 10 head of the file.
entrezID   dbXrefs
399518   VGNC:VGNC:43937|Ensembl:ENSCAFG00845006432
399530   VGNC:VGNC:40330|Ensembl:ENSCAFG00845002136
399544   VGNC:VGNC:48329|Ensembl:ENSCAFG00845029798
399545   VGNC:VGNC:41624|Ensembl:ENSCAFG00845011460
399653   VGNC:VGNC:44358|Ensembl:ENSCAFG00845001610
403152   -
403157   VGNC:VGNC:38980|Ensembl:ENSCAFG00845013158
403168   VGNC:VGNC:37699|Ensembl:ENSCAFG00845014982
403169   VGNC:VGNC:54270

My input gtf file
#!genome-build ROS_Cfam_1.0
#!genome-version ROS_Cfam_1.0
#!genome-date 2020-09
#!genome-build-accession GCA_014441545.1
#!genebuild-last-updated 2020-10
X       ensembl gene    24550462        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; gene_source "ensembl"; gene_biotype "protein_coding";
X       ensembl transcript      24550462        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; tag "Ensembl_canonical";
X       ensembl exon    24552206        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; exon_id "ENSCAFE00845128634"; exon_version "1"; tag "Ensembl_canonical";
X       ensembl CDS     24552206        24552226        .       -       0       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; protein_id "ENSCAFP00845021332"; protein_version "1"; tag "Ensembl_canonical";
X       ensembl start_codon     24552224        24552226        .       -       0       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; tag "Ensembl_canonical";

So the final goal is to add entrezID to the modified gtf file keeping the format same as it was in the input file.
Here in my small gtf input the gene_id is only for ENSCAFG00845015183 as due to space constraint i didn't put more lines.
So now how do I map my reference file entrezID to the gtf file gene_ID as this pattern ENSCAFG which is common in both the file.
Now in my mapping file I dont have entrezID for ENSCAFG00845015183 since due to first 10 row i have used to show.
So how do I map these two file where the modified gtf fill will have new attributes as entrezID

Please [edit] your question and i) add an example of your input file (you linked us to a directory, show us the actual file) and ii) show what output you expect from that example. We need to know what ID you need (gene? transcript? protein? something else?) and the exact format you are trying to achieve.

okay let me add more details to it

@terdononstrike i have updated my question please do have a look and let me know how do i start

@terdononstrike i found this tool https://agat.readthedocs.io/en/latest/tools/agat_sq_add_attributes_from_tsv.html which can do what i need but if i have to do it in shell script how to go about it?


text:
Statistical methods suitable for DE analysis of non coding RNAs
I am currently working on DE analysis of coding as well as non-coding features of A. thaliana using the edgeR package.
Is the negative binomial method that is normally used for the DE analysis of genes suitable for non-coding features such as miRNA, lncRNA, snRNA, snoRNA, pseudogenes, transposons, uORFs and 5' and 3' UTRs? Or are there special statistical methods that are more suitable for modelling each of these, like I have read that the zero-inflated model used by lncDIFF can better model lncRNA counts? I see that my miRNA counts matrix is abundant in zero counts and very low counts of genes. Should I use lncDIFF for that?
Also, among the various DE test methods in edgeR such as exactTest(), glmFit() and glmLRT(), glmQLFit() and glmQLFTest(), which is/are the one(s) most suitable for running DE tests on such non-coding features?

There is no question that lncDIFF is the approach to use, especially if the alternative is performing edgeR directly, and thats demonstrable (below). I don't think anyone would question the legitimacy of your approach. If you look at Fig 1 for the paper describing the test, here there is no question about it edgeR used directly does not perform for this task at all. The paper describing lncDIFF is here,
https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-019-5926-4
You present lncDIFF as an option within edgeR. This I don't know, but lncDIFF  can be directly downloaded from GitHub, https://github.com/qianli10000/lncDIFF.

The one thing you have to be careful about is you are using plants, whilst these tests are designed around humans, see below for details.

How it's tested: GLM methods versus an Exact test, e.g. Fisher's. Both are all valid tests. However, GLM methods are preferred  over an exact test, in my statistical experience.
First, you need to obtain a fit before you can perform a test, which ever GLM you choose.
x <- model.matrix(...)
y <- estimateDisp(...)
fit <- glmFit(y, x, robust=TRUE)
lrt <- glmLRT(fit, contrast=c(-1,1,0,0,0))

The basic issue is GLM likelihood ratio test versus quasi likelihood test, which the right approach? The clear answer is QLT because the authors of lncDIFF developed the test around QLT (see the paper). That shouldn't be general advice outside lncDIFF
The correct approach is to try both, if there major differences thats an issue: I would say in that scenario is consider  QLT - but only for lncDIFF, as described. If there are no major differences you are good to go.

Thanks. I read the paper.  That answers the first part of my question. What about the last part? Also, among the various DE test methods in edgeR such as exactTest(), glmFit() and glmLRT(), glmQLFit() and glmQLFTest(), which is/are the one(s) most suitable for running DE tests on the other non-coding features?

Hi @ArkajyotiBanerjee upvotes and/or "accepted answers" are much more preferred to thanks: you can upvote. My advice is to open (a) separate question(s) for the remaining issues because this would clarify what is being sought. Multi-part questions are generally discouraged.


text:
Filtering criteria for non-coding features with very low counts
I am trying to do DE analysis of non-coding features of A. thaliana. I find in the miRNA and lncRNA counts file that they are abundant in zero counts, and most of the non-zero counts are very low. Now, for the DE analysis of coding features, an inbuilt filtering method to remove the extremely low counts is available in edgeR, called filterByExpr(). Is this function suitable for only those features with high counts like coding sequences, or is it suitable for non-coding features with low counts and low library sizes too?
If not, then what kind of filtering criteria should I use for filtering out genes with insignificant counts? Or should I filter out only the ones with row sums 0?

Partially answered before. I am assuming this is a lncDIFF calculation. For the lncDIFF paper it is unlikely to work as the following figure displays,
https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-019-5926-4/figures/4
However, don't take my word for this, test it. You should test the option explicitly via glmQLFTest and glmLRT. If you perform either of these tests with or without the filter the better fit (one not rejected) is the correct answer. I strongly predict the filtered data will be rejected and unfiltered data  (hopefully) will pass.
The question of what is better glmQLFTest versus glmLRT is discussed in the previous answer, here Statistical methods suitable for DE analysis of non coding RNAs
Just keep in mind this is plant data and these tests are developed for humans.


text:
Gene Ontology (GO) analysis and classification of bacterial whole genome annotations
I am new to Gene Ontology (GO) analysis and would like to classify and visualize the annotations. Here is the background and the steps I have taken so far:

I annotated the bacterial whole genome using Prokka. I performed

Interproscan on the annotated genome using the Galaxy server.

The Interproscan results are in a tabular format, and column 14
contains the GO IDs.

Example GO IDs from column 14:
column 14
GO:0003677
GO:0003677
GO:0005524|GO:0016887
GO:0005524|GO:0016020|GO:0055085
GO:0003723|GO:0046872
GO:0004521|GO:0004534|GO:0006396|GO:0008270|GO:0090501
etc

I would like to classify these annotations into the Cellular Component, Biological Process, and Molecular Function categories, and then plot a graph similar to the provided image.

I am comfortable using either R or Python for this analysis. Could you please provide guidance on how to achieve this using the preferred language?
Thank you in advance for your help!

`library(clusterProfiler)
library(org.Hs.eg.db)
data <- read.delim("A3_InterProScan.tabular", header = TRUE, sep = "\t")
go_ids <- data[,14]
gene_symbols <- mapIds(org.Hs.eg.db, keys = go_ids, keytype = "GO", column = "SYMBOL")
gene_symbols <- gene_symbols[!is.na(gene_symbols)]
enrich_result <- enrichGO(gene = gene_symbols, OrgDb = org.Hs.eg.db)`

Hmmm... well I know how I'd solve this and its a lot of code. I'd download the GO database (I assume it's downloadable) and screen the names from there. Then there's a pandas based solution to obtain frequencies and plot it. That strategy isn't a single question. Again, its possible there's something published.

@M__, Thank you, If other easy way is ther thats fine too, Kindly provide me a code, I would be grateful to you, Thanks in advance

The site is not a coding service. As the question stands published solutions are good answers. There's one possible solution, but I've never really used it especially in this context of the question, so can't recommend it.

A python version of this question has been posted here: https://stackoverflow.com/questions/76614627/

Thanks @RamRS provision of code is always good. In my opinion this is four separate questions in terms of code 1) getting the gene function from GO; 2) obtain a frequency of each gene function; 3) arranging gene function into groups; 4) plotting the resulting frequency.


text:
Algorithm for counting neighbors
I need an algorithm to process protein chains. The volume size would probably not exceed 30x30x30 angstrom^3.
Say, I have a point cloud of 50,000 points in a 3D space.
I want to count the number of neighbors each point {(x1,y1,z1), (x2,y2,z2), (x3,y3,z3), ... ...,(xN,yN,zN)} has at radii {r1, r2, r3, ..., rM}.
Which algorithm should I use?
Additional information:

Do you need to do this count online/sequentially or all at once?
All once.

What restrictions might there be on the radii?
The radii will be 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0 angstrom.

Do you need to do this just once in this particular case or do you need a general algorithm?
I will use this to process protein sequences. Not DNA/RNA, though.

For a general algorithm, what's the likely range of the number of points?
The range would be 50,000 at max (although, the maximum length of available protein is 38,000 residues).

Are you expecting there to be relatively few or relatively many neighbors?
I am expecting relatively few neighbors, in the range of 0-50.

NOTE: I prefer a hash-based algorithm as it doesn't require additional data structures like KD tree, or octree.

The usual method is dividing the space into cells, in your case 7Å cubes, assigning atoms to cells, and searching for neighbors only in 3x3x3=27 cells. This method is described in Wikipedia as cell lists, although that name is not as widely used as the method itself.
If you analyze experimental protein structures note that they may have alternative conformations (alternative locations of the same atom), and that the structure in the file represents an asymmetric unit.

K Dimensional Tree (KD-Tree) data structure might help:
https://www.biorxiv.org/content/10.1101/605907v2.full.pdf

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
Assigning multiple GO IDs to an ID: Should I split or use the first or last ID in the hierarchy?
I have a dataset that contains IDs and associated Gene Ontology (GO) IDs. Some of the IDs have multiple GO IDs assigned to them, and I'm unsure how to handle this situation. I would like to determine whether Should I split all Go id or I should use the first or last GO ID in the hierarchy when assigning the GO ID to the ID.
Here is an example of my data:
3667 FEHLGCFF_05007                                             GO:0006637|GO:0047617
3673 FEHLGCFF_05919                                  GO:0016020|GO:0022857|GO:0055085
3684 FEHLGCFF_00369                                             GO:0022857|GO:0055085
3695 FEHLGCFF_00369                                  GO:0016020|GO:0022857|GO:0055085
3696 FEHLGCFF_00369                                                        GO:0022857
3702 FEHLGCFF_06641                                                        GO:0016627
3707 FEHLGCFF_06641                                                        GO:0003995
3721 FEHLGCFF_04414                                                        GO:0046872
3725 FEHLGCFF_02968                                  GO:0003700|GO:0006352|GO:0006355

How should I handle cases where an ID has multiple GO IDs assigned to it? Should I split all Go id or Should I use the first GO ID, the last GO ID, or is there another approach I should consider?
Any guidance or best practices would be greatly appreciated. Thank you!
Here is My code using R, it it Ok
library(tidyverse)

# Read the go-basic.obo file
go_data <- readLines("go-basic.obo")

# Extract the GO ID and Gene ID from the unique_data dataframe
unique_data <- unique_data %>%
  separate_rows(GO_ID, sep = "\\|") %>%
  select(V1, GO_ID) # Replace with the correct column names from your unique_data dataframe

# Create an empty dataframe to store the results
result <- data.frame(Gene_ID = character(), GO_ID = character(), Name = character(), Namespace = character(), stringsAsFactors = FALSE)

# Process each line in the go-basic.obo file
for (i in 1:(length(go_data) - 1)) {
  line <- go_data[i]

  if (grepl("^id:", line)) {
    go_id <- strsplit(line, " ")[[1]][2]

    # Check if the current GO ID exists in the unique_data dataframe
    if (go_id %in% unique_data$GO_ID) {
      # Filter the unique_data dataframe for the current GO ID
      filtered_data <- unique_data %>% filter(GO_ID == go_id)

      # Extract the name and namespace information
      name <- ""
      namespace <- ""

      while (!grepl("^is_a:", go_data[i])) {
        if (grepl("^name:", go_data[i])) {
          name <- strsplit(go_data[i], "name: ")[[1]][2]
        } else if (grepl("^namespace:", go_data[i])) {
          namespace <- strsplit(go_data[i], "namespace: ")[[1]][2]
        }
        i <- i + 1
      }

      # Add the information to the result dataframe
      temp_df <- data.frame(Gene_ID = filtered_data$V1, GO_ID = go_id, Name = name, Namespace = namespace, stringsAsFactors = FALSE) # Replace with the correct column names
      result <- bind_rows(result, temp_df)
    }
  }
}

# Print the result dataframe
result
nrow(result)

What do you want to do with this data? GO terms of deep down the hierarchy also imply the upper nodes (genes annotated in the olfactory sensors term are also part of the biological process term).

@ I I have seen in papers, they classify in GO terms and give good bar graphs with frequency, whith three terms, but I got this kind of confusing, so I am confused which one is correct, for publishing paper

There are three branches of GO terms, BP, biological process, CC, cellular components and MF molecular functions. All three are correct, for publishing, it only depends on what do you want to focus on. They are independent of each other, that's why they are split/classified in three terms.

I think I should delete this question

Dear @umar I am afraid that deleting a question with an upvoted answer is not permitted. It is fair to point out both llrs and RicardoA.'s answer are reasonable contributions. Also RicardoA's comment about how you obtain GO assignments is good.

You might want to look at Revigo to simplify multiple GO terms into their parent terms. I've used rrvgo extensively to simplify enrichment results - this is the R implementation of Revigo.
It's not clear how exactly you are obtaining these GO assignments. If these are the significant hits for each feature, I'd take the union of all significant GO terms (i.e. everything in your second column combined) and pass it to rrvgo. Then, you can select the parent term associated with each row. I don't guarantee each row will be collapsed to a single term, but it should make whatever you are doing downstream much more tractable.


text:
is it possible to count cell in Violin plot in seurat?

I tried to quantify the number of cells of Fzd9+ Macrophage.
how can I know the cell count?
thank you

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

Does @haci's response answer your question? If so please mark it as "accepted". The question appears to contain very little detail and is easy for a misunderstanding to occur in these situations.

This can be a good starting point: table(your_seurat_object@meta.data$the_column_that_has_the_labels_of_interest)
Basically you will be counting the number of rows (each single cell is a row) in the meta.data slot of your Seurat object for each level of your column of interest (for example cell type annotations).


text:
Feasible to find genetic variations of two samples using RNAseq data?
I have bulk RNAseq data from two strains of mice from Jackson Lab: C57BL/6 and B6.SJL.
The former expresses a Ptprc-b allele and the latter expresses a Ptprc-a allele.
They are supposedly isogenic at all other loci. But to the best of my knowledge, there is no empirical evidence of this.
My question is, how would I go about finding all the homozygous single nucleotide variations (if there are any) in the coding sequences of the two strains from my RNAseq data.
I should mention that the RNAseq data are paired-end illumina reads ~150bp long.

Maybe try MMseq or MMseqs2?


text:
DE analysis tool and method for non-coding features
I am currently working with the non-coding features of A. thaliana, and trying to get the DE features. Among the three DE test methods in edgeR such as exactTest(), glmFit() and glmLRT(), glmQLFit() and glmQLFTest(), which is/are the one(s) most suitable for running DE tests on miRNA, snRNA, snoRNA, pseuogenes, transposons and 5' and 3'-UTRs?
Should I be using any other DE analysis tool instead of edgeR for some specific non-coding feature, like there's lncDIFF for lncRNA?

Regarding what are the differences between these tests, this answer by one of the authors might be helpful: https://support.bioconductor.org/p/76790/

The post here is useful ... https://support.bioconductor.org/p/76790/
ExactTest is really cool and completely different from glmQLFit glmFit etc ... you'll need two controls to check it works (it could just reject everything). Its described here https://rdrr.io/bioc/edgeR/man/exactTest.html
There is a concern, some of the lncDIFF data looked more like a lognormal to me in the publication, but I didn't read it in sufficient detail. Anyway your data will probably approximate to a negative binomial distribution and it's statistically tricky to  assess whether two negative binomial distributions are significantly different ... the ExactTest will do this.
Is the total output of the control significantly different to the test sample? ... the ExactTest will tell you. You want p < 0.05, preferably a lot lower than 0.05 to be convinced. However, you need to compare two control samples against each other, where p > 0.05 for the subsequent tests to be valid.

So really nice BUT you need to check whether it is a negative binomial FIRST. Thats what the other tests are doing. Please do check whether lncDIFF identified a negative binomial, because if it ain't ExactTest will not work, they could be using a different distribution.
The other tests you mentioned are different to the ExactTest because they simply assess whether the data is negative binomial. The distribution is the frequency of each RNA transcript BTW. edgeR's lnfDIFF requires that distribution and thats how it can assign probabilities to the difference between the frequency distribution of two identical transcripts under difference experimental conditions (or other conditions). So you need a test to say "yeah I'm a negative binomial" (assuming it is negative binomial), then everything else via lncDIFF is then valid.
First you've got to fit you data to a model ...
y <-glmGLMFit(...)

... is the function and that uses general linear models to fit the negative binomial. The fit will result in something called "a residual" and that is very important, the smaller the better. Testing whether that residual is small enough for a negative binomial to be declared is down to ...
glmLRT(y ...)

This is a likelihood ratio test, to test the observed likelihood against the model (p > 0.05). You can also use
gof(y ...)

this is goodness of fit, as the link says the residual follows a chi-squared distribution, i.e. its tested on the chi2 tables and thats easy to understand. Again you need a probability > 0.05 to declare a fit.
HOWEVER, the lnfDIFF guys developed a quasi-likelihood to assess the distribution
z<-glmQLFit(...)

Fits the data via quasi-likelihood, I honestly have no empirical idea how that works.
glmQLTest(z ...)

tests the fit to the distribution. You need a probability > 0.05. Again I have not empirical idea whats going on here.
ANY non-significant result is good, doesn't matter if its GLM with gof, GLM with LRT, GLM-QL (whatever that is) because all of these are solid tests.  Which would I use? The guys who developed lncDIFF used QL (quasi-likelihood) to assess the fit (might be a good idea). However, statistically I honestly don't know whats going on, so I would also look at the GLM and fit it against the goodness of fit - because thats easy to understand.
What do I do to assess the fit? I use machine learning, but it's way more complicated than anything here.

From the comments.
It's important to note the context. The OP is performing non-coding RNA work, lncDIFF is much better than normal edgeR for this type of data described here Statistical methods suitable for DE analysis of non coding RNAs.
lncDIFF is now available for R here
... so the methods of edgeR should be available to it. glmQLFit() is a general function in R. Providing the lncDIFF object will be accepted by the tests, by performing everything in R this should solve this.
I would advise using lnfDIFF's QL directly, i.e.
y<-ZIQMLfit(...) # fit the distribution
ZIQML.LRT(y ...) # test the fit

The manual is here
https://github.com/qianli10000/lncDIFF/blob/master/lncDIFF-manual.pdf

I hope that makes sense these are pretty advanced stats so it can be tricky to understand them. Whilst I can teach code, I don't teach complex stats outside machine learning.

So... is edgeR a good tool for analysing the DE of noncoding features? I mean since lncRNA has a separate specialised tool for its analysis, which they claim is better than edgeR, is edgeR fine if I want to analyse noncoding features?


text:
Issue with Differential Gene expression analysis with Deseq2
I am using the Deseq2 package to perform a DE analysis on multiple factors.
So, basically to perform the analysis I consider patients characterized by the Brugada Syndrome type 1 and other patients that instead are controls.
In order to have equal numbers of PZ and CTRs I included in my analysis all the PZ of type 1, that in this case are 8 samples, and I selected 8 random CTRs among the full dataset the CTRs that I have for this batch (in total are 40).
To select the random CTRs I used the sample_n function:
random_CTR_samples_IGATECH<-sample_n(samples_IGATECH_ALL_CTR,8)

Anyway I created the design matrix that looks like this including also the Age variable as well as the Gender one since I want to test also the effect of those in my data, (NON_INDUCIBLE)= CTR:
       Sample BrS_BASELINE_PATTERN Gender Age
1  DG1340                    1      M  44
2  DG1383                    1      F  35
3  DG1453                    1      F  55
4  DG1455                    1      F  53
5  DG1473                    1      M  33
6  DG1557                    1      M  45
7  DG1578                    1      F  57
8  DG1590                    1      M  42
9  DG1380       NON_INDUCIBILE      M  20
10 DG1381       NON_INDUCIBILE      F  52
11 DG1382       NON_INDUCIBILE      M  49
12 DG1565       NON_INDUCIBILE      M  30
13 DG1594       NON_INDUCIBILE      F  42
14 DG1620       NON_INDUCIBILE      M  40
15 DG1633       NON_INDUCIBILE      F  62
16 DG1652       NON_INDUCIBILE      M  21

and I filtered the Count matrix for the samples present in the design one. So I have
head(Counts_IGATECH_ALL_1_random_CTR)
            DG1340 DG1383 DG1453 DG1455 DG1473 DG1557 DG1578 DG1590 DG1380 DG1381 DG1382 DG1565 DG1594 DG1620 DG1633 DG1652
ENSG00000114771      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0
ENSG00000129673      0      1      0      0      1      0      0      1      0      1      0      0      0      0      0      0
ENSG00000144452      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0
ENSG00000143921      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0
ENSG00000173357      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0
ENSG00000161103      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0

When I run the DESeqDataSetFromMatrix function I get this error:
converting counts to integer mode
factor levels were dropped which had no samples
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

  Please read the vignette section 'Model matrix not full rank':

  vignette('DESeq2')

and the function i ran is:
dds_IGATECH_1_random_CTR<- DESeqDataSetFromMatrix(countData = round(Counts_IGATECH_ALL_1_random_CTR),
                                               colData = design_IGATECH_1_random_CTR,
                                               design = ~ Age + Gender + BrS_BASELINE_PATTERN)

I revised my design matrix multiple times and tried to change also the combinations of the CTRs to have a different one and see if that would solve the problem but it didn't work, and I cannot figure out what it can be.
Thank you in advance for any suggestion.
Update
Design matrix:
dput(design_IGATECH_1_random_CTR)

structure(list(Sample = c("DG1340", "DG1383", "DG1453", "DG1455",
"DG1473", "DG1557", "DG1578", "DG1590", "DG1376", "DG1379", "DG1573",
"DG1575", "DG1579", "DG1598", "DG1627", "DG1633"), BrS_BASELINE_PATTERN = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), levels = c("1",
"2", "NON_INDUCIBILE"), class = "factor"), Gender = structure(c(2L,
1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L), levels = c("F",
"M"), class = "factor"), Age = structure(c(23L, 15L, 32L, 31L,
13L, 24L, 34L, 21L, 29L, 27L, 30L, 10L, 19L, 28L, 18L, 37L), levels = c("17",
"20", "21", "22", "23", "24", "25", "26", "28", "29", "30", "32",
"33", "34", "35", "37", "38", "39", "40", "41", "42", "43", "44",
"45", "46", "48", "49", "50", "51", "52", "53", "55", "56", "57",
"58", "61", "62", "63", "64", "69"), class = "factor")), row.names = c(NA,
-16L), class = "data.frame")

The NON_INDUCIBLE data are different since i have generated a new design matrix again, but the format is the same

This typically happens when Age is a factor rather than numeric column, so each unique row is nested with itself. Design ~BrS_BASELINE_PATTERN + Gender is fine, there is no nesting, so this narrows it down to the Age column. Cannot check as you do not provide reproducible data via dput().

Thank you for your fast response, by converting it to a numeric column it works. Still have this note though

the design formula contains one or more numeric variables with integer values, specifying a model with increasing fold change for higher values. did you mean for this to be a factor? if so, first convert   this variable to a factor using the factor() function the design formula contains one or more numeric variables that have mean or standard deviation larger than 5 (an arbitrary threshold to trigger this message). Including numeric variables with large mean can induce collinearity with the intercept. Users should center and scale numeric variables in the design to improve GLM convergence.

This now goes into the realm of how you should design your analysis which is not what the question is about. There are many previous threads via google, e.g. on Bioconductor support site on numeric covariates. Numeric covariates like age assumes that there is a relationship between age and expression, which is questionable. It can be that old is different than young, but a 40-yo is unlikely to have double or have expression as compared to a 20-yo, but this is what this would assume with a numeric covariates. Scaling age is an option, or categorizing (young-middle-old).


text:
Why is there gene without annotated location?
I came across this gene Pepc (peptidase C)
NCBI gene database has a page for it here https://www.ncbi.nlm.nih.gov/gene/109616
But it does not contain any information on its coding sequence or genomic location. It just says it sits on chromosome 1.
Can someone familiar with how the reference genome is annotated explain why this is?

I agree that this is an odd gene: no transcript, but it has mouse strains and supporting literature.

I like the question a lot. You can't have this sitting on arguably (sure EBI would disagree) the worlds most comprehensive molecular database, not as comprehensive as it thinks it is. This is one for escalation with NCBI and the authors. If this was EBML Emily would step in. Definitely upvoted, I would recommend site members upvote to help the OP recover their bounty costs. I am aware this sort of thing goes on ... I've encountered in a completely different area.

Let me know if you intend to escalate it, because I will otherwise.

@M__ Thanks for the comment. I'm not sure what you mean by escalating. Do you mean reaching out directly to NCBI?

Yes, you can approach the authors but it never results in success (I've tried - not on this example).

This gene/protein appears to have been described exclusively as an isozyme / electrophoretic mobility assay of a protein.
See for example this PepC mutant page:

The Pepc locus controls electrophoretic mobility of the enzyme peptidase C of red blood cells. This LT/J congenic strain carries the Pepca allele from C57BL/6 that determines a slow anodally migrating electrophoretic band.

The relevant publications describing mutants are from the 70s and 80s, e.g. here, when there were no genomic resources for mice.
However, linkage mapping of electrophoretic variants could be used to establish linkage to loci on chromosome 1, apparently.
Very likely PepC is encoded by some annotated gene on chr1, but no one knows which one exactly because they don't really know much about the protein other than its electrophoretic variability. But that's just a guess by me.
When the databases were set up, they had to accommodate 1) mouse mutant strains, and 2) published articles referencing the protein, so they did their best.
However, no one ever followed up to link everything together, so we just have these orphan entries.


text:
How evaluate biclustering algorithm over real datasets containing labels?
I have a set of real data containing labels. However, in some biclustering algorithms (e.g., Cheng and Church Algorithm), originally authors applied gene expression data like Yeast without having labels. My question is that what are appropriate evaluation measures in this scenario?

Welcome to the site. I understand machine learning (ML) but the question is not clear. The key Issue is the use of links. We will not click two  external links to understand the question, so it needs a summary of the intended algorithm is needed within the question. Biclustering sounds to me to be unsupervised learning and that would explain the difference. More generally, explaining what data you have and what the labels are is usually helpful.

Biclustering or co-clustering is an unsupervised learning task with which both rows and columns in a data matrix can be placed in concurrent clusters whereas non-overlapping clusters are generated in clustering algorithms.
Data set - a resource of regulatory elements (REMs) contains features like Gene ID, REMID, Genomic location, Predicted function (acts as a label having values such as activation, repression), Model score (normalized absolute binary logarithm of the P-value).

Briefly evaluation measures, which I presume means "features", are determined using supervised learning using a technique called "feature selection". Unsupervised learning generally is used to identify categories for a supervised learning algorithm. The orthodoxy has changed.
"Unsupervised learning feature selection"

To identify suitable features using unsupervised learning a priori then genome location which must have a major outcome, thereafter it would be a case of deleting a given feature and re-evaluating the clusters. If the clusters disintegrate then that is a key feature (not very cool but would work for datasets with limited numbers of features).
Weighted unsupervised learning is popular, this is basically "feature selection", i.e. which features contribute to the clustering. This is basically the formal version of point 1. Personally I have had little success with this approach, i.e. the results were not informative, but it does work. The results are not discriminatory between different conditions, which is what I was seeking.

Regarding labels, labels are not used in unsupervised learning approach - except weighted unsupervised learning - hence their presence or absence usually makes little difference. The only "label" used is the sample type, i.e. the "id" and not the "feature". Without intrinsically understanding the experiment it is difficult to forward an accurate a priori assessment.