text:
How do I plot clusters based on DNA sequence alignment?
I have a group of about 1050 sequences that have been aligned using blastn. I am now trying to create a cluster plot of the sequences in order to identify similar sequences, something like this:

However, I am having trouble finding a tool to do this. My professor recommended using ClustalW, and so I have spent plenty of time searching for a tool that can take a sequence alignment and plot them such as in the linked image. However, I have not been successful in finding what I need and am beginning to think I am missing something or something is going over my head. If anyone could assist with finding the tool I need I would very much appreciate it. Thanks so much!

I think one of the reasons you struggle is that clustering DNA sequences is not a clearly defined task. In general intention of clustering is to reconstruct, or approximate the relatidness of the DNA seqeunces.
If all these sequences are homologous you can do a multiple sequence alignment (using MAFFT or Clustal) and build a phylogenetic tree (using something like IQ tree). If your prof suggested ClustalW, this is probably what they meant. Evolutionary speaking, that's the optimal method for getting relations of sequences, but they must be "alignable". Also with this many sequences, it will be computationally quite heavy.
For within species relatedness, phylogenetic reconstruction is sometimes impossible, because recombination breaks down the linkage of individual alleles. However, if all of the sequences are long enough, you can do a principal component analysis on SNPs, for example using something like plink. The plot you showed actually looks a lot like PCA (it even has the explained variability on axes).
Finally, regardless of their relationship, you can always do a dendrogram based on a similarity matrix, but for that, you need to define similarity. That will be usually something kmer-based. I am not following this field very closely lately, but there are some recent papers about kmer-based clusterings, like this one.
When looking for the paper I also run into this one, which proposes genome clustering using "Inter-nucleotide Covariance". No idea how useful/sound it is. Just found the abstract quite interesting.
Long story short, it really depends on what you are after and what are your sequences! If you can build a reliable tree, it is the most powerful way to determine the relationships among sequences.


text:
Modeling number of reads mapped to a gene
I am looking for a probability distribution of a number of reads mapped to a particular gene in metagenomic sequencing (NGS, shotgun, likely illumina).
Naively one could model it via a binomial (or Poisson) distribution as:
$$
P_n = {N\choose n}\pi^n(1-\pi)^{N-n},
$$
where

$n$ is the number of reads mapped to the gene
$N$ is the sequencing depth (total number of reads)
$\pi\propto l$ is the probability of read mapping to the particular gene, proportioal to the gene length (possibly adjusted for the read length)

This is however complicated by such intricacies of sequencing as:

finite length of the fragments, which can be longer or shorter than the gene length
paired-end sequencing, which means one has to be careful when counting evens where only one or both reads map to a gene
something else?

Is there an established model for this case? Could you recommend reading on the subject?

This is quite a similar question to this one, but I notice that you're specifically asking about a probability distribution rather than differential expression. That sounds like gene dispersion and/or shot noise; a good reference for me on that has been the documentation for DESeq2:
http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
The gene dispersion issue is discussed in more detail in their 2014 paper:
https://doi.org/10.1186/s13059-014-0550-8

We model read counts K ij as following a negative binomial distribution (sometimes also called a gamma-Poisson distribution) with mean μ ij and dispersion α i . The mean is taken as a quantity q ij, proportional to the concentration of cDNA fragments from the gene in the sample, scaled by a normalization factor s ij , i.e., μ ij =s ij q ij . For many applications, the same constant s j can be used for all genes in a sample, which then accounts for differences in sequencing depth between samples. To estimate these size factors, the DESeq2 package offers the median-of-ratios method already used in DESeq. However, it can be advantageous to calculate gene-specific normalization factors s ij to account for further sources of technical biases such as differing dependence on GC content, gene length or the like, using published methods, and these can be supplied instead.

source: Stack Exchange

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

sort: votes


text:
Convert .smp downloaded from NCBI CDD into HMM(.hmm) file
Is there a method or tool to convert a .smp file from CDD (Protein Domain File) into a HMM file? I downloaded the CDD files from the FTP site and am more familiar with domain searching with the HMMer suite.

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


text:
What are the best practices for comparing .seg copy number outputs from iChorCNA to VisCapCancer?
I have shallow whole genome sequence data that I used the tool iChorCNA to find the mean copy number segmentation per region. The tool produces a seg file and can be seen here: https://github.com/broadinstitute/ichorCNA
I have also run these same samples to higher coverage on a targeted panel and used the tool VisCap Cancer to find the mean copy number segmentation. The tool can be found here and also produces a seg file: https://github.com/pughlab/VisCapCancer
I am looking to compare the outputs of each tool to see how close they overlap with each other. The problem is that iChorCNA produces outputs in bins that span the genome per million bases, and VisCap Cancer is only looking at the sites for targetted probes. These probes are anywhere from 1.3 million basepairs to 260 million basepairs in length.
What I am looking to do is segment the VisCap output into bins that match the bins used by iChorCNA. I want to then take the sum of the mean value given by iChorCNA at each bin that overlaps with the VisCapCancer results, and divide it by the total number of bins that fall in the range of the results to find the average iChorCNA reading per VisCap reading. Then, I can do a t-test, ANOVA or other measure to see how closely these two numbers are to each other.
I was wondering if anyone has directly compared copy number results from two data sources such as this before, such as shallow whole genome sequencing against a targeted panel, and if they have any advice as to how to go about this most effectively - including if it would be better to do the comparison using R, python or another language, and best practices for setting up the comparison most effectively.
Thank you for your time


text:
What is the meaning of split read?
I want to use rna seq data to later perform functional tests on fusion genes.
so before that I need to filter the "best results" (of rnaseq) for deciding which candidates I actually want to try in the lab.
one criterium is the split read. (among coverage, discordant mate and confidence)
In some cases the split read (or discordant mate/coverage) of one of the two genes that fusioned is "0".
Does that automatically mean that the gene is no longer a candidate for functional testing?
Unfortunately I only got the results in an excel data sheet and don't know which tool was used to gain data.

For discordant mapped reads, a position of 0 for the pair of a mapped read usually means that the pair couldn't be mapped anywhere. This would be expected for fusion events because some reads could originate from across the fusion site (i.e. not being fully contained in either "unfused" gene), leading to a sequence that doesn't match anything in the reference.


text:
Parse RNA variant effect annotations ("r." format)
I've got annotations for splicing variants in a format like this (this is one variant):

Variant: NM_004092.3:c.88+5G>A
Effect: Retention; r.87_88ins1_88+10:p.(Ala31Glufs*23)

I want to extract which acceptor/donor sites are lost or gained (in this case lost causing intron retention). In other words, the genomic location of the sequence that is retained in mRNA.
I parsed the variant HGVS ID using ensembl VEP (like this) without problems, giving me the genomic location of this variant. However, VEP seems to not like "r." annotations for the effect: I tried ensembl VEP again with RNA, protein, and both annotations, none of which worked.
Is there a tool that can parse these annotations for me? I've got too many to do it by hand.

Note on RNA splicing
RNA splicing and acceptor and donor sites:
An introns contains: a 5' donor site, a branch site located near the 3' site and a 3' acceptor site. Both the donor and acceptor  are required for splicing along with the branch site.


text:
How to decide between standardized mean difference (SMD) and log2fold change (log2FC) before doing meta-analysis?
I want to do meta-analysis of alpha diversity between test and control groups across 10 different populations. I want to see what is the summed effect (pool effect) of intervention between the two groups for alpha diversity. For that, first I need the intervention effect on test group's alpha diversity for individual population dataset. How should I calculate this? Should I calculate this by difference in the log transformed mean values of the groups (log2FC) or should I calculate SMD for each study?
Thanks

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

what is the measure of alpha-diversity? counts or proportion?


text:
How to detect bedtools version used by pybedtools (in order to correctly preserve strand information when merging gtf records)
I have issues working with pybedtools due to strand information being lost after extracting and merging transcript coordinates from a gtf file.
It seems that the solution to preserve strand information depends on the version of bedtools, so I would like to have my code detect bedtools version and adapt its behaviour accordingly.
On a machine with bedtools v2.27.0:
ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).saveas()
print(my_genes)

I       ensembl transcript      4119    10230   .       -       .       gene_biotype "protein_coding";gene_id "WBGene00022277";gene_name "homt-1";gene_source "ensembl";gene_version "1";p_id "P18791";transcript_biotype "protein_coding";transcript_id "Y74C9A.3";transcript_name "Y74C9A.3";transcript_source "ensembl";transcript_version "1";tss_id "TSS21161";
X       ensembl transcript      17714968        17718714        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8395";transcript_biotype "protein_coding";transcript_id "cTel55X.1b";transcript_name "cTel55X.1b";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";
X       ensembl transcript      17714968        17718717        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8874";transcript_biotype "protein_coding";transcript_id "cTel55X.1a";transcript_name "cTel55X.1a";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";

merged = my_genes.sort().merge(s=True).saveas()
print(merged)

I       4118    10230
X       17714967        17718717

(Missing strand information)
merged = my_genes.sort().merge(s=True, c=7, o="distinct").saveas()
print(merged)

I       4118    10230   -
X       17714967        17718717        +

On a machine with bedtools v2.26.0:
ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).saveas()
print(my_genes)

I       ensembl transcript      4119    10230   .       -       .       gene_biotype "protein_coding";gene_id "WBGene00022277";gene_name "homt-1";gene_source "ensembl";gene_version "1";p_id "P18791";transcript_biotype "protein_coding";transcript_id "Y74C9A.3";transcript_name "Y74C9A.3";transcript_
source "ensembl";transcript_version "1";tss_id "TSS21161";
X       ensembl transcript      17714968        17718714        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8395";transcript_biotype "protein_coding";transcript_id "cTel55X.1b";transcript_name "c
Tel55X.1b";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";
X       ensembl transcript      17714968        17718717        .       +       .       gene_biotype "protein_coding";gene_id "WBGene00007068";gene_name "cTel55X.1";gene_source "ensembl";gene_version "1";p_id "P8874";transcript_biotype "protein_coding";transcript_id "cTel55X.1a";transcript_name "c
Tel55X.1a";transcript_source "ensembl";transcript_version "1";tss_id "TSS18767";

merged = my_genes.sort().merge(s=True).saveas()
print(merged)

I       4118    10230   -
X       17714967        17718717        +

merged = my_genes.sort().merge(s=True, c=7, o="distinct").saveas()
print(merged)

I       4118    10230   -       -
X       17714967        17718717        +       +

(Duplicate strand information)
Given that I work on linux systems, how can I detect bedtools version from within python?
Is there a better way to do what I want?
(Ideally, I would like to have my resulting bed consist in chrom, start, stop, name, score, strand fields, where the score would just be "." and the name would be the content of itvl["gene_id"])

Edit (05/05/2021): Regarding my last comment above, I found out that it is easy to convert gff into bed using pybedtools.featurefuncs.gff2bed (tested with bedtools 2.30):
from pybedtools import BedTool
from pybedtools.featurefuncs import gff2bed

ids = ['WBGene00022277', 'WBGene00007068']

def in_ids(itvl):
    return itvl[2] == "transcript" and itvl["gene_id"] in ids

all_genes = BedTool("genes.gtf")
my_genes = all_genes.filter(in_ids).each(gff2bed, "gene_id").saveas()
print(my_genes)

I       4118    10230   WBGene00022277  .       -
X       17714967        17718714        WBGene00007068  .       +
X       17714967        17718717        WBGene00007068  .       +

merged = my_genes.sort().merge(s=True, c=[4, 5, 6], o="distinct").saveas()
print(merged)

I       4118    10230   WBGene00022277  .       -
X       17714967        17718717        WBGene00007068  .       +

This ensures proper bed format.
However, with bedtools 2.26, s=True makes strand information appear as a 4th column, so it is more complicated to obtain correct column order after merging.

pybedtools assumes that bedtools is in your path and bedtools itself will return the version with bedtools --version. So:
import subprocess
subprocess.check_output(['bedtools', '--version'], text=True).strip().split()[1]

For me that returns 'v2.30.0'.


text:
Finding targets genes of putative transcription factor after homer motif enrichment
Steps I have done

I have two cell types to compare one is stem cell and leukemic stem cell
Two sets of data for each condition RNA seq and ATAC seq.
Did differential expression as well as differential accessibility analysis.
From the differential accessibility analysis I'm getting the regions which I used for motif analysis using homer.
The result of de novo motif enrichment

Question
So this results shows possible transcription factor binding motif. Now how do i find the target genes of this putative transcription factor. Do i need to see the literature or what other ways to do it. As of now I don;t have chip seq data. So is it possible to use expression and accessibility data to find out the possible targets of this putative transcription factors, which in turn can give me lead for my downstream validation.
Any suggestion or help would be really appreciated


text:
How is the odds ratio of disease risk conferred by a 1-standard deviation increase in PRS calculated?
One standard deviation from the mean is commonly used to calculate a polygenic risk score for GWAS, e.g. human genetic disease.
Why is this a common metric, for example why not 2-SD or 1.96 SD as in the normal distribution and what is the formula?
Here's a figure taken from this paper (10.1186/s13073-020-00742-5) to help guide the explanation:

The "OR/SD" metric is somewhat arbitrary. It's just convenient that it is 1 SD; there is no exact reason. If you wanted to make it OR per 2 SDs, you could just multiply the OR by 2.
Imagine you have a distribution from 1st to 100th percentile of PRS, and you have a linear association between percentile of PRS and probability of disease. How would you quantify this? OR/SD is a quick and easily interpretable way to say how sensitive the PRS is. In other words, how "fast" the OR increases as you increase PRS.
The method for obtaining this metric usually involves converting your raw scores into Z-scores. The slope of your regression where Y is OR and X is Z-score I think would correspond to OR/SD.


text:
Calling for one specific SNP from multiple sequencing runs
Just learning the basics of bioinformatics and bash shell scripting.
Is there a correct way to search for a specific allele with one or even multiple sequencing runs? I am searching for prevalence of the C2416T allele in the SARS-CoV-2 genome within the ORF1a RNA polymerase gene as a mini project.
I downloaded a bunch of SARS-COV-2 sequencing runs, I think I need to convert them to bam file type than I am not sure how to call for one specific SNP, especially in a loop for 5-7 runs. Can this be done with bash?
I was looking into variant calling from a bam file type using vcf and I can't find a way to call for one specific snp on multiple bam files at once.
Ideally I would want to filter the variant out of the .bcf files or use grep since I sort of understand those the most and I’m not working with a large amount of data or a human genome.

Getting SNPs from whole-genome data is usually done via mapping sequencing reads on the reference (bwa-mem or bowtie are two popular mappers) and then using one of variant calling tools (like freebayes, samtools, or GATK). However, I would also look up if there is something virus-specific.
According to this answer, you could get variant calls of a specific region using samtools if you give it a bed file that contains the region of interest (three columns chr    start    stop format). For example using human reference
samtools mpileup -uDl regions.bed -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - > RAL_samtools.raw.bcf

You could also specify a bunch of files, not just two and it would do the variant calling on them all.

Answer from @atpoint [transferred from comment]:
Alternatively, you can use samtools view to only retain reads that overlap the region you are interested in and then feed these into the mpileup. IIRC there was something with that -l option in mpileup that caused trouble but I do not remember what it was. I think it still called all variants based on all reads and then simply filtered for the target regions in -l. The view approach would probably be much faster. Just keep that in mind as an alternative.


text:
Use Htslib to extract specific region of vcf file
I am working in c++ using the htslib library to parse different vcf files as well as fasta files.
My goal is to read in a specific region from a fastafile and then extract all the VCF sites within that region from a seperate vcf file, in this case only containing chr14. I can extract the fasta region in several ways, one of which is to use the function faidx_fetch_seq() and define a region of chr14:20000000-22000000.
Then for the bcf files i am able to iterate through the entire file and extracting the different information regarding the reference, alternative allele and so on. With a simple example shown below.
void faVCF(const char* fastafile,FILE *fp){
  htsFile *test_bcf = NULL;
  bcf_hdr_t *test_header = NULL;

  test_bcf = bcf_open("chr14.vcf", "r");
  test_header = bcf_hdr_read(test_bcf);
  bcf_hdr_set_samples(test_header, "20000000-22000000", 0);
  bcf1_t *test_record = bcf_init();
  while(bcf_read(test_bcf, test_header, test_record) == 0) {
              std::cout << bcf_hdr_id2name(test_header, test_record->rid) << "\t" <<
                          test_record->pos << "\t" <<
                          test_record->n_allele << "\t" <<
                          std::endl;
          }
}

However, i doesn't work to extract all the variant sites in a specific region. I tried bcf_hdr_set_samples with the argument "20000000-22000000", to extract that given region from my chromosome. As far as I understand this should be able to subsample both samples and region queries using the header information, and then upon initiation extract that region, however this doesnt work.
So my question is: Can anyone help me with what function I can use to extract all the vcf sites from the region chr14:20000000-22000000.?

For iterating over reads in a region in a single BCF file, you want to use the iterator interface: bcf_itr_queryi() or bcf_itr_querys(), together with bcf_itr_next(). For a VCF file, it appears you need to use the similar TBX iterator interface: tbx_itr_queryi() or tbx_itr_querys(), together with tbx_itr_next(). See HTSlib's tabix.c for reasonably straightforward examples of both.
Alternatively you could use the routines in htslib/synced_bcf_reader.h, which encapsulates both and is what bcftools mostly uses.


text:
Bioconductor, genefilter() returns NULL dimensions, is my filtering process wrong?
I'm using R and bioconductor in order to conduct some gene analysis on an Affymetrix dataset (GSE173360) but I'm having trouble at filtering genes. I'm trying to create a dataset called small.eset by filtering my dataset but small.eset returns NULL dimensions. I was wondering if I'm doing some obvious mistake?
The dataset I'm using is GSE173360 and I'm sharing my condensed code (without plots and quality checks).
# Library and Data Imports ------------------------------------------------
library(affy)
library(affyPLM)
library(simpleaffy)
library(limma)
library(RColorBrewer)
library(genefilter)
setwd("/Users/Documents/R Scripts/biostat_geo/data")
Data <- ReadAffy() #read all CEL files in working directory

# log2 transformation -----------------------------------------------------

ex <- exprs(Data)
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0)
if (LogC) { ex[which(ex <= 0)] <- NaN
exprs(Data) <- log2(ex) }

#background correction

eset <- expresso(Data, bgcorrect.method = "mas",
             normalize.method = "loess",
             pmcorrect.method = "pmonly",
             summary.method = "mas")

# filtering

esetunlog <- 2^exprs(eset)
ffun <- filterfun(pOverA(0.20,10))
t.fil <- genefilter(esetunlog,ffun)
small.eset <- (log2(esetunlog[t.fil,]))


text:
Finding 2000 bp sequences that occur more than once in the genome
I am looking for a software that can identify long sequences (1500-2000 bp) that occur more than once in my genome of interest. A k-mer counter like KMC could have worked, but KMC has a max limit of 256 bp. If I can club contiguous 256-mers from KMC, that would work too, but I don't know how to do that. Any ideas?

MUMmer includes a repeat-match tool which can do exactly that:
$ repeat-match
USAGE:  repeat-match  [options]  <genome-file>

Find all maximal exact matches in <genome-file>

Options:
 -E    Use exhaustive (slow) search to find matches
 -f    Forward strand only, don't use reverse complement
 -n #  Set minimum exact match length to #
 -t    Only output tandem repeats
 -V #  Set level of verbose (debugging) printing to #

Note that a "maximal exact match" is just a maximal repeat.

This looks ideal. Will try this out and update. Thank you!

Okay I tried it out and it worked, so thank you! Ended up using mummer instead of repeat-match so as to compare between (and not within) sequences


text:
Convert genbank file into tbl file
I want to submit some phage sequences to NCBI. So I want to transform a genbank file into a tbl file. Can someone help me please?
Thanks

Is tbl a Sequin file format equivalent ??  https://chlorobox.mpimp-golm.mpg.de/GenBank2Sequin.html

It's a perfectly valid question


text:
DESeq2 multiple treatments, multiple time points, multiple cell lines
Note: this question has also been asked on Bioconductor Support
I know there are a lot of questions asking similar things here on this forum and I checked the vingette and did a lot of other research but I still have a hard time wrapping my head around it.
In my position, data is unfortunately just put in front of me and I'm asked for pretty pictures. Even making it clear to the wetlab people in the lab that we need replicates is an issue.
I have multiple cell lines, multiple time points and multiple treatments:
Cell lines: CL1, CL2, CL3 Time points: 6h, 24h Treatments: T1, T2, T3, Control
For each cell line and each time point, there are 3 different treatments plus a control. 3 replicates for each sample -> 72 samples
What I want is actually quite simple; I want to measure control vs each treatment at each time point in each cell line. No testing across time points or across cell lines.
My first thought was to separate the data into 6 different data sets (CL1, 6h | CL1, 24h | CL2, 6h | CL2 24h | CL3, 6h | CL3, 24h) and simply do the DE analysis separately, but I read that this is not the way to go. In the end it would also be nice to get a normalized count matrix with all data normalized together, for PCAs and similar.
I hope I made it understandable. How do I design the DESeqDataSet?

As Mike Love alluded to on the Bioconductor pages, it's best to have specialised help for datasets that are this complex; doing analysis on anything which has multiple compounding factors is not simple. It's also a good idea to talk through the analysis before doing sequencing, so you can discuss things like having six replicates per sample (including controls).
In any case, what I did for something similar to this (eventually, after lots of discussion with biologists about possible models) was to combine all the non-tested factor columns into a single factor. Something like this:
meta.combined.df <- meta.df %>%
  mutate(lineTime = paste0(cellLine, "_", timePoint))

And then subset the data based on this statistic:
counts.sub.mat <-
   counts.sub.mat[,meta.df$lineTime == desired.lineTime]
meta.sub.df <-
   subset(meta.combined.df, lineTime == desired.lineTime)

Then use ~ Treatment as the design formula (assuming "Treatment" is a field containing either the treatment type or "control").
But... this is explicitly treating each line/time point as a separate experiment, which may not be what you want. There are complications associated with batch effects and outliers that need to be talked through, because you could end up identifying "differential expression" that is not reproducible.

When comparing one subset of samples to another subset of samples, you generally do not have to make a subset dds object.  It's preferable to leave all the samples in for dispersion estimation.  (Though if you wanted to look for linear trends across two treatments and multiple time points, you will have to either omit the unwanted treatments, or do some trickery with the model matrix).

Feel free to add this information to your answer. In my case, subsetting worked best for the data we were working with.

To compare a subset of samples to another subset of samples, make one column of column data which will distinguish them all the way you want (so probably join cell line, treatment, timepoint), and make that column your design, and then when you call results, you specify the pair you want to contrast, one pair at a time.


text:
Problems about value log2(IP/Input) less than zero in ChIP-seq?
I have a question about the normalization for ChIP-seq. I used CPM to normalize my bam files of each IP and Input. Then I calculate the coverage of gene bodies for all genes on the genome. I have WT and different mutants. Take H3K4me1 as an example, I first calculate the matrix over gene bodies for H3K4me1 IP and Input. Then I used the matrix of H3K4me1 divided by the matrix of input. then log2(H3K4me1/Input) and plot this value over entire gene bodies. I found all the regions are lower than 0 and it seems that my IP does not work. In fact, this is not the case since I have checked the antibody's efficiency via ChIP-qPCR and it works well. Moreover, the distribution of H3K4me1 is consistent with the published data. But the problem is the profile of log2(H3K4me1/Input) is lower than 0.  I don't know how to deal with such problems? I would appreciate it if someone can give me some suggestions on that.  Thanks!

Seconding what Devon says:

Use a proper normalization technique such as TMM from edgeR. For applying this to bigwig files you could follow my tutorial at biostars.org.

My two cents beyond that:

I would not use enrichment over input as input is hard to normalize and can be sparse depending on sequencing depth. Therefore, you might get biased enrichments that do not well reflect reality, especially if counts in IP or input (of both) are low. I personally only use input for peak calling, and for plots like this then compare IPs among each other.

So the strategy would be to first normalize your bigwigs (I guess you use these for the plots?) using the strategy as linked above, and then use these robust CPMs to make the plots without any input samples, only using the IPs. You can use MA-plots to explore whether samples are properly normalized. I made a little code suggestion here in this biostars post on how to do this starting basically from a matrix of raw counts using edgeR. This is basically (from the TMM normalization standpoint) the same as in the first link, but then doing pairwise MA-plots to see whether TMM properly centers the majority of points along y-axis being somewhat zero. Normalization is key in any experiment, therefore I always recommend to use these kinds of plots for diagnostics.

Dear ATpoint, sorry for replying to your message so late and thank you very much for your suggestions!  Yes, you are right, I first want to normalize the raw bigwig files once I get the normalized bigwig files, I will try to calcualte the matrix of genebodies for genes interested in each sample.  I will have a try to use TMM method you recommended to see whether it works on my data or not!  I will let you know once I got the results.   Thanks!

Hi, ATpoint, I read the link you sent to me. I have a question about the regions in the matrix . If I understand properly,  for example, I have 5 strains and each strain have a peak region list, you means that pooling the peak coordinates of each sample as one region list,  then use the feature count to calculate the coverage of each region for all samples, if I am wrong, please correct me. I am looking forward to your reply! Best.

This is a common issue caused by the normalization method not being particularly robust (most likely it's using the mean signal for normalization, which will be a bit upward skewed when you have regions of high enrichment in the ChIP). The basic ways to fix this are:

Use a more robust normalization (e.g., the median signal, or even just manually select un-enriched regions to calculate the normalization factors)
Manually tweak the normalization factors so the background (i.e., log2 ratio in regions where you don't expect any enrichment) is 1. That will generally make the regions you're plotting more biologically realistic.

Note that there can be a negative enrichment in some areas, due to things like exclusion of histones, so don't be overly zealous if you manually tweak the normalization factors.

Dear Devon,  Thanks for your kind reply and valuable suggestions! I have tried to use the median value of the non-enriched region to calculate the scaled factor of each IP vs their input. I found it works for all modifications but failed for H3 IP. I guess this is due to the peaks for H3 are widely spread on the genome and when I use MACS2 to do the peak calling for H3, it won't work efficiently. Do you have any suggestions for H3 IP normalization? I am looking forward to your reply.  Thanks!

You may have to manually tweak normalization values to get reasonable results.

What I did: first did peak annotation, then extracted the genes(name it as gene list 1)  with enriched peaks. Then I calculated the matrix of the gene body of genes (not included in gene list 1 ), further calculate the median value for each sample (for both IP and input).  In the peak annotation step, I tried Homer and Rgmatch, unfortunately, the results of these two software are quite different.

Therefore, I recalculated the median value given the coverage of the non-peak region (just escaped the peak annotation step). It seems the results of using Homer were much closer to the results directly given the scaled factor for the non-peaks region.  I checked the default parameter of Homer, they defined TSS as from -1kb to +100bp and TES as -100bp to +1kb. Since S.cerevisiae. genome is not big like human or mouse, I am not sure how much effect using the default parameter.   I am pretty interested in the 500bp up-and downstream of gene bodies.

The Homer source code is based on perl and I haven't used Perl before.  Do you have any other suggestions? Thanks！


text:
Alignment tools on desktop for peptide alignment sequence and how to make annotations on alignments look presentable?
I have 10 transcripts where I need to do 5 peptide sequence alignments (2 transcripts per alignment). I need to do this in order to indicate which amino acids lead to which domains.
My questions are two:

Is there a tool I can use from my laptop that aligns peptide sequences? Because results from the clustal omega peptide alignments web service only last two weeks, they're not ideal for replicability in a paper.

The ultimate aim is to use those alignments and domain classifications to a project I am writing and I am trying to make it as presentable and scientific as possible. Any suggestions on how to make peptide alignments and annotations on peptide alignments (which indicate protein domains) look presentable for a paper or if a paper has already shown peptide alignments in a presentable way that would be highly appreciated.

Could you specify what do you mean by peptide, i.e. are we talking about oligopeptides, 5-20 AA, commonly referred to as simply peptides? Or do you mean protein chains. For the former, domain makes no sense, while for the latter, many tools better than Clustal Omega exist (as Omega is for large datasets, most common being Muscle). About annotations on protein, the term features is used and there are a few feature viewers... So you may want to specify webapp, Python or R here.

Hi Matteo, thank you for your response.
I want to differentiate between classical class I and non classical class I MHC molecules in a model organism using well conserved structural features within classical MHC I molecules (eg intradomain disulfide bridges within class I a1 domain). To do this I need to align the peptide sequences of transcripts of that model organism with transcripts of classical MHC I human genes and see which genes from those transcripts have those structural features. So I am looking for a tool to align and annotate those structural features in a replicable manner

Additionally if you have any other bioinformatic ideas /computational analyses on how to find classical MHC Class I molecules on an organism which is poorly annotated that would highly appreciated!


text:
Shannon's Entropy and DNA
I have been working with Shannon's entropy (SE) and DNA lately, and referring to the formula and concept of SE, I just wondered whether one should in the case of DNA use each nucleotide as having a probability to show up as 0,25 or if it was better to use the probability based on the frequency in the genome (let's say the human genome).
In the case of the latter being the better, is it interesting to have those frequencies at the scale of the genome, the chromosome, the sub segment ?
I'm interested in calculating the entropy relative to a short window (the length of an Illumina read, about 125-150bases), using the proportion of each base on a given sequence length. The idea is "how much information is given" in this sequence, and by information I mean how much different bases do I have.
It might that I'm missing some basic mathematics knowledge here, don't hesitate to point it out so I can understand it better
cheers,

I think that genome-wide probabilities (i.e. not 0.25 for each base) should be fine. For a given segment of DNA, the chance that the sequence doesn't appear anywhere else in the genome is fairly high.
There is lots of patchy non-randomness (e.g. at centromeres, STRs, and other tandem repeat regions), but there's not much point in modelling that because it's not random.

Thanks for your answer. So if I understand it well, not taking equivalent probabilities but the human genome bases probability is a better way ?

@Bratten Yes; I've made a slight update to my answer to clarify this.


text:
Base abbreviations other than ACTGU in sequence file
Short question:
where can I find what F and J and < and - and 7 mean in a Sequence file?
Longer version:
I was looking at this sequence file I got from NCBI (it's an sra file? if that is the correct terminology) and I used the sratoolkit to print out a few spots
looking at them I expected to see ACGU (this was meant to be an RNA sequence), but what I see are like

AAFFFJJJFJJJJJJJJJJJJJJJJJJJJJFJJF<JJ77AJJJJ-JJJJJFF7FA--F7FAJJ7<AJ-FJA<J<JJFA<<77F-JJJJJJF<JF<7-FF-A-----A-77JJ7FJJF-FFFAJ<<-FJJFAJAFFAFFFFF<7A7JFJJ7AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJFJJJJJJJJJJJJAA<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJAAFFJJFJJJJJJJJJJJJJJFJ-FJF

Where can I find a reference to help me "understand" this?
I mean I assume - means missed or unknown? or is it meant to be already aligned and these are deletions?
I appreciate any pointer to get started. as you must have observed I am a total noob.
I got the file from: https://www.ncbi.nlm.nih.gov/sra/SRX8156420[accn]
Then I did a
 ./fastq-dump -X 5 -Z ERR4557575

Follow-up question:

Where do you look for the quality scores ranges and such, I guess as @winnie2k suggested I should convert this to FASTQ and make my life easier, but still wanted to understand the file better.

Regarding the follow-up question: SRA files are a red herring. Nobody I know cares about the SRA file format ¯\(ツ)/¯

I usually convert SRA files to FASTQ before doing any further processing, but the tooling for converting SRA files is cumbersome.
To simplify things, I recommend you download the FASTQ file directly. The ENA currently has a nicer interface for doing these things. Here is a link to download the FASTQ files for the SRA sequencing run associated with the SRA experiment you are describing: https://www.ebi.ac.uk/ena/browser/view/SRR11588968
A description of the FASTQ format can be found here: https://en.wikipedia.org/wiki/FASTQ_format


text:
10bp sliding window quality score
I am having trouble creating a 10 basepair sliding window that goes through seq and trim them all 10 basepairs that have an average quality score (qual) less than 15. I am trying to approach this by reading in the file and seperating the sequences and quality scores, then read in the 10 basepair window in the sequence.
#!/usr/bin/envpython
with open('fastaQ.fq') as f:
    lines = [line.rstrip() for line in f]
    seq= []
    qual= []
    tenFrame= []

    #Seperates the sequences and the quality score of the file
    for i in range(len(lines)):
        if lines[i].startswith('@'):
            seq.append(lines[i+1])
            qual.append(lines[i+3])
            i=i+4

    #
    for i in range(0,len(seq)):
        tenFrame.append(seq[i][i:i+10])
        if len(tenFrame)%10 == 0:
            print(tenFrame)```

I'm not sure what you mean by "trim". The code below will discard any reads that have a 10bp window with a mean phred score below your cutoff. I'll ammend the post to also include a solution that "trims" basepairs on both ends until the window achieves an average score of above the cutoff. It would be ideal if you could provide a minimal reproducible example that shows what your code does and what you'd like it to do.
As far as the code goes, I'd first suggest using a library like Biopython to parse the Fastq file. It will make your life a lot easier. For your code, it looks like seq actually contains all the sequences, so len(seq) is the number of reads in the file. However, when you iterate over i in range(len(seq)), it looks like you're using i to both index which read to iterate over as well as get the 10 base pair window. You need an additional index variable to obtain the window:
for window_start in range(len(seq[i]) - window_len) + 1)
    window = seq[i][window_start:window_start + window_len]

Additionally, some pythonic tips for the future:

If a range starts from 0, don't specify it as an arugment i.e. use range(x) as opposed to range(0, x)
Try to avoid anything along the lines of range(len(x)). In this case, you could just iterate over the sequences directly: for s in seq: or, if you wanted the index, you could enumerate the sequences via for seq_idx, s in enumerate(seq)
Use good variable names. seq implies a single sequence, yet is a list of sequences. Same with qual

Here is some code that will retain only sequences that have all windows with a passing mean phred score
from Bio import SeqIO

read_file = "/file/here.fq"
w = 10
cutoff = 15
valid_seqs = []
with open(fq_file) as fq_fd:
    seqs = SeqIO.parse(fq_fd, "fastq")
    for seq in seqs:
        scores = seq.letter_annotations["phred_quality"]
        mean_scores = (np.mean(scores[i:i+w]) for i in range(len(scores) - w))
        if any(mean_score < cutoff for mean_score in mean_scores):
            print(f"{seq.id} does not pass the filter!")
            continue
        valid_seqs.append(seq)


text:
Cellranger gives error
I am trying to run cellranger but I get fastq permission denied error
fi1d18@RBGO-Server2:/opt$ cellranger-atac count --id test --reference refdata-cellranger-arc-GRCh38-2020-A-2.0.0 --fastqs /data/fi1d18/Downloads/FASTQ_scATAC/ --chemistry ARC-v1
error: couldn't open MRO file
    Caused by: Permission denied (os error 13)
fi1d18@RBGO-Server2:/opt$

I have added the path by
export PATH=/opt/cellranger-atac-2.0.0:$PATH

And by vim I have also added the path
fi1d18@RBGO-Server2:/opt$ vim ~/.bashrc

Any help please?

Does your user have permission to read inside opt ? mro are configurations files used and read by cellranger during each steps. If they cannot be access, then the pipeline should fail.

Two of those command line parameters should be accessible locations (see here). Could you please add the output of the following commands:
ls -l /data/fi1d18/Downloads/FASTQ_scATAC/ | head
and
ls -l refdata-cellranger-arc-GRCh38-2020-A-2.0.0 | head
I suspect that the refdata link should be an absolute path, rather than a relative path, i.e. something like this:
cellranger-atac count --id test --reference /opt/refdata-cellranger-arc-GRCh38-2020-A-2.0.0 --fastqs /data/fi1d18/Downloads/FASTQ_scATAC/ --chemistry ARC-v1


text:
Computing Maximum Agreement Subtree of two unrooted phylogenetic trees with PhyloNet?
I am using PhyloNet to compute the MAST of two unrooted binary phylogenetic trees but I am not getting the output that I should. I think the issue might be with the way I am writing the trees in the input, but  I don't know.
For example: if I input T1=((((1,5),6), (2,4)),3) and T2=((1,2), (((3,5),6),4)), the MAST returned is ((1,2),3), whereas if I draw the trees, I get the MAST ((3,6),(4,2)).
I think this might happen because somehow it assumes a root, but I have no idea how to fix it or if that is really the issue. Any guidelines for computing MASTs of unrooted trees are appreciated.


text:
Getting exact mismatch data from multiple sequence alignment
TL;DR Can anyone point me in the direction of a tool that can be used to input different viral substrain reference sequences (fasta) as well as sequence data (fastq/bam/vcf) and find the most likely substrain of the sequenced data?
I am doing multiple sequence alignment of different viral subtype references and would like to get the exact positions and base each time there is a mismatch. I have tried using different tools. For example Clustal Omega can almost give me what I seek, but I have trouble getting the exact position and nucleotide for the mismatches (the full thing is ~8000bp long):
ACTACAATAATCCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAATTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG

*********** *********** ****** *****************************

E.g I would like something like:
Pos 12:
C
T
T
T
T
T
T
T
T
T
T

For each position with a mismatch.
A problem is that the sequences are slightly different lengths so something like a variant call format on the "main" strain might work, but I'm not sure how I can do this on multiple files of this length.
The end goal is actually to be able to call which substrain is present from NGS data. I imagine using variant calling on the "Main" strain, but getting from there to calling which substrain is in the data is a bit of a long way for me. If someone have any suggestions on tools for this, please enlighten me :)
Thank you very much!

I found a solution.
For variant calling on multple sequence alignment I can recommend snp-sites.
For anyone facing a similar problem of wanting to find substrains, here is an overview of useful tools:
"Computational Methods for Strain-Level Microbial Detection in Colony and Metagenome Sequencing Data"
(https://www.frontiersin.org/articles/10.3389/fmicb.2020.01925/full)

I have partly answered this for know! For anyone interested i use kalign to align all 10 strains to a fasta format, then I use snp-sites to output a vcf file that can show each mismatch in accordance to the first strain in the fasta. New problem is that the alignment "stretches" the length of this first strain, because some of the other strains are a bit longer. This means that the position of variants are with reference to the longest strain of the 10 and not the strain I want. It can work for now, but any solution to this, would be helpful.


text:
Biological datasets for polythetic classification
Question: What are some instances of polythetic datasets in biology? In particular, I am looking for a dataset to benchmark a machine learning algorithm optimized via episodic training.

On polythetic classification
The conventional definition of a conceptual or monothetic class is that its members must possess certain properties in common [1]. For example, in the MNIST database of handwritten digits, images can be classified into 10 distinct classes (i.e. digits) where all instances of any given class possess the same set of attributes (e.g. zeros are typically written as ellipses).
In contrast, a polythetic class encompasses members that might share a number of characteristics, but none of these attributes is essential for membership of that group or class. For example, one could group green zeros together with red ones in one polythetic class, and red zeros together with green ones in another one, as shown below.

Are there any instances of biological datasets with similar properties? I am particularly thinking of genetics data (e.g. single nucleotide polymorphisms) where perhaps a certain mutation A counteracts the effect of mutation B? Any pointers are welcome.

Hi @rvinas, this question seems far too general to be useful. There are, for example, plenty of examples of gene sets that have negatively correlated expression under different conditions, and heritable syndromes that have different physiological appearances depending on the parent they were inherited from. Could you please be more specific about what sort of data you're looking for.

I may not be understanding exactly what you are talking about, but it seems to me that nearly any set of classes of interest must be polythetic. Most groups of comparable things differ in some characteristics but also have characteristics in common, and none of those characters is 100% identical with membership in the class.
Some examples off the top of my head:

Biological species / taxonomic groups (classify according to phenotypic characters)
Biological species / taxonomic groups (classify according to molecular sequence e.g. DNA)
Cell types (classify according to cell epitopes)
Cell types (classify according to transcriptional profile)
Ecosystems (classify according to species profile)
Disease state (classify according to any one of many molecular profiles)
Members of population groups (classify according to genetic variation)

etc.


text:
Is it possible to correct bond order and conectivity problems in an sdf file according to mol2 files?
I was doing a Virtual Screening Search using Pharmer software. The .sdf output file contains a column in mol format and another column with Root Mean Square Deviations. I can see some bond order and connectivity problems in the mol column. I do have a mol2 file that contains the original molecules and I was wondering if I can use these mol2 files, to fix this. Can rdkit work in this case?

RDKit has a one line command to fix the bond order (docs)
 mol = Chem.MolFromMolFile('wonky.mol')
 template = Chem.MolFromMol2File('template.mol2')
 newMol = AllChem.AssignBondOrdersFromTemplate(template, mol)

The problem is that the two need to be identical. I see a carboxylate in your screenshots, so you most likely did the correct thing and docked with the appropriate protonation for pH 7, i.e. as the conjugate base with a negative formal charge. If one does not then RDKit will fail.
Also, RDKit is "not great" with mol2 files. If you have the original SMILES that is better. If the problem is atom names, using a PDB is a better call. If you want to assign atom names, this blog post of mine may help.


text:
How can I speed up FASTA ssearch output?
I'm trying to predict targets for miRNA sequences, and so far after reading some papers for other tools that do this, I've found that I should be using FASTA's ssearch36 for this.
I'm starting with a small file of 500 miRNA sequences, searching a reference DB of 50k sequences. When using ssearch, I get the following output at the very bottom after ~85 minutes of running:
Total Scan time: 294.000 Total Display time: 4681.000
I am piping the output to a file, and the result file is ~13 Gb. It is taking this long even when I pipe the output to >nul so I don't think it's my disk speed.
I'm wondering if there's any way to speed this up, since it seems very unreasonable to me that the actual search itself only takes ~5 minutes but yet taking over an hour to simply generate output. Is it a bug or am I interpreting it incorrectly?

are you using the parallel or MPI version?  If so, you can add more threads with the -T option.  However, this probably won't help since the default is to use as many processors as are available


text:
How to make a video of a gromacs gmx dump standard output with VMD
I did gmx trjconv on GROMACS trr file to generate a .xtc file, and then used gmx dump and saved its standard output in a file. Then I tried to make a movie with VMD using the file. As a result, it does not work properly. The error in the terminal is as follows.
gromacsplugin) Cannot read header fromm '/Users/kondouakira/mount/Lab/1ns_sim/md_0_1_noPBC3.trj', no error

I use following command.
gmx trjconv -s md.tpr -f md.trr -o md_0.trr -pbc mol -ur compact -b 1 -force
gmx dump -f md_0.trr > md_0.trj

No need to use gmx dump, you can simply load a xtc or trr trajectory in VMD using:
vmd path/to/gro_file.gro other/path/to/trajectory_file.xtc

As always you can also use the GUI to first open your gro file, then right click on the molecule and select Load data into molecule before selecting your xtc or trr file.
After putting the representations as you see fit, if you want to export this as a movie (GIF, MP4, others...) you can go in Extensions -> Visualization -> Movie Maker. See Trajectories and Movie Making VMD official tutorial for more details.

Thank you for your answering! But I want to draw a video using a trj file.There is also a trj file option in the VMD file format. Therefore, I would like to know how to perform this to create a video.

I don't think that the output of gmx dump can be read by VMD: the goal of this command is to provide a human-readable output, not to convert one trajectory format into another.
trj files are NOT the output of gmx dump. trj files were used in "old" versions of GROMACS to store trajectories (see trj file format documentation of GROMACS 5), and this is what VMD is awaiting, not the output of gmx dump. This is also why it tells you the header is wrong.
In recent versions of GROMACS trj files are no longer supported, probably as they were not in a binary format, hence took too much disk space (see Trajectory files formats in GROMACS 2021 documentation). Converting a .trr in a .trj was done by gmx trjconv, but I don't think this is possible anymore (see gmx trjconv supported trajectory formats in GROMACS 5 vs GROMACS 2021).
Why would you absolutely need gmx dump output to be read by VMD? The only reason I can see is a handmade script doing modifications to your trajectory. In this case I would recommend to instead look into some molecular dynamics python libraries, like MDAnalysis, PyTraj or MDTraj, which will also allow you to make minute modifications to the trajectory file while still reading it in its original format. You can also look harder into GROMACS tools, as they are quite complete.


text:
How can I associate centrality score of genes(vertices) in a GRN with the probabilistic model in terms of cause and effect?
I'm learning to build Bayesian Network Based Modeling of Gene Regulatory Network of Cancerous Cellular Cells. So far, I have learned to form a Bayesian Network based on nature of genes(i.e Functional Genes and TFs). However, I cannot figure out how to give probability to each node/vertex. Is it possibile to somehow iterpret the centrality score of each vertex an associate the meaning with the probabilistic model in terms of cause and effect?
or can you atleast refer me a relevant research paper?
I believe I would be worth mentioning that I'm fairly new to this.
Here's my data of TFs and Sources.

TF
Target
Model of Regulation
References

ABL1
BCL2
Repression
11753601

AHR
FOS
Unknown
10478842

APC
BIRC5
Unknown
11751382

AR
ERBB2
Activation
21613411

AR
ERBB2
Unknown
21741601

Following is the sample of my data for network analysis, mostly containing centrality scores some graph properties.

SUID
name
BetweennessCentrality
ClosenessCentrality
ClusteringCoefficient
Eccentricity
EdgeCount
Indegree
IsSingleNode
NeighborhoodConnectivity
Outdegree
PartnerOfMultiEdgedNodePairs
selected
SelfLoops
shared name
Stress

7228
TP53
0.06611272
0.46456693
0.04141865
4
75
50
FALSE
14.40625
25
11
FALSE
0
TP53
39572

6788
AR
0.03605041
0.40410959
0.03164983
4
60
46
FALSE
13.38181818
14
5
FALSE
0
AR
20498

7029
MYC
0.03342125
0.35014837
0.03280053
5
81
74
FALSE
12.93589744
7
3
FALSE
0
MYC
19140


text:
Alignment to verify plasmid from forward and reverse Sanger sequencing reads
I am analyzing Sanger sequencing data to verify whether or not the correct protein coding sequence was inserted during cloning.
For each sample, I have both forward and reverse Sanger sequencing reads that cover the coding sequence.
Would it be better if I first align the forward and reverse reads to get a consensus, and then align the consensus to the reference?
Or should I just align each read directly to the reference?

Align all reads from a sample at the same time IMO. This will give you all of the information at once, and it should not be difficult. It should be very obvious whenever there's a disagreement between reads or between reads and reference.
I would also recommend inspecting chromatograms directly at all sites of interest (such disagreements, insertion sites, possible errors or mutations).


text:
Imputation of dog genotype
I am an undergrad working on a set of genotype matrix of dogs, it's in csv file format with information on only the each sample as rows, chromosome and its position as columns, and the rest of the matrix is filled with 0,1,2 and missing values.
So I am wondering what would be the best method to impute the genotype with all the data I get? I looked up Impute2, but it seems to require the reference data(recombination map) which I don't have...
Really appreciate whoever answers my question!!

I have never inputted genotype data, but I am pretty sure a recombination map for dogs exists. Check for example here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5100850/

I've done imputation on dog data before - so I can hopefully help. However, I think we need some more information from you. How many individuals and SNPs do you have to begin with? Are the dogs the same breed or mixed?

Impute2 is a great, if dated choice. Another great choice is beagle. The advantage of beagle in this case is that it is not as strict about requiring a recombination map. Without a map, results will be poorer, but perhaps good enough.
Do make sure that your sample size including a reference panel is at least 200 in any case.


text:
Annotatation of transposable elements / transposons / SINE sequences in human
I need to annotate human genome with different types of information, amongst them: transposable elements.
I read about repeatmasker, but what I understood is that it estimates the probability of a sequence to be a transposon based on the sequence and not based on curated database.
My university does not have access to RepBase (GIRI). I was hoping finding some possibility with biomart but it seems not possible.
What would be a way to access curated annotation of all kinds of transposable elements (SINE, Alu seq..) within a given genome/genome coordinates?
My goal is to estimate the link between some measurements we have made on the genome (binned) ( and Shannon entropy is one of them, and probably the less interesting regarding TE). Eventually, I'd like to be able to estimate that measurement X has a correlation with the number of Alu repeats present in the neighbourhood.
What would be the best way to achieve this?

repeatmasker is meant to annotate/mask a genome given a library of known transposable elements. From the project description, I assume you have not done any sequencing on your own. Then, instead of annotating, I think you should find already computed annotations of TEs.
The alternative to subscription-based repbase is dfam. And Just looking at their webpage they seem to have an annotation of all individual TEs in hg38. I bet my shoes there will be plenty of TE annotations in human in genome databases too.

I'll check it out, thanks. Indeed, for now we have sequencing data, but given our analyses we're willing to use them "as is" as they have also been already aligned. Thanks !

Thanks @kamil-s-jaron it was definitely helpful. I managed to intersect with what I wanted. My only issue is that I don't really get which type of TE is it. There's plenty (>500) family names in my subset and it's not spontaneously possible to say that BC200 is a type of alu repeat for instance. As you seem to know a lot about this topic, would you have a recommandation to link my findings to a type of TE? (dfam proposes a classification which could actually be helpful, but I couldn't find a way to call back this information ?) thanks

Honestly, I don't know. But it's also a completely new question, would you mind opening a new one? Someone else might be able to help you.

Thanks, will do that


text:
Differentiating molecules based on peptide sequence? How to annotate?
I want to differentiate between classical class I and non classical class I MHC molecules in a model organism using well conserved structural features within classical MHC I molecules (eg intradomain disulfide bridges within class I a1 domain). To do this I need to align the peptide sequences of transcripts of that model organism with transcripts of classical MHC I human genes and see which genes from those transcripts have those structural features. So I am looking for a tool to align and annotate those structural features in a replicable manner
Additionally if you have any other bioinformatic ideas /computational analyses on how to find classical MHC Class I molecules on an organism which is poorly annotated that would highly appreciated!

MAFFT is able to incorporate structural information in alignments: https://mafft.cbrc.jp/alignment/software/dash.html

Please consider the answer below and upvoting/accepting where appropriate.

Knowing almost nothing about the biology, I think that the most straightforward way to do this is:

train profile HMMs on your different classes of protein sequences (in human maybe, but ideally including MHCs from a large number of organisms)
translate your transcripts into proteins and run them against your trained profile HMMs
confirm that the profiles adequately discriminate the two classes and adjust if needed.

For an example of how to do this (on MHC), you can see here.
These approaches are then easily extensible to further organisms, as you have the profile HMMs.
It's possible that a simpler approach would consist of just finding hits against human MHC and then just checking for the existence of sequence patterns corresponding to those structural motifs, but I doubt that would work well in any significantly divergent organism.


text:
How to create this function for lapply?
I have two lists of dataframes. List x and list my_list. Each having four dataframes and their names are same for the two lists, i.e., bray-curtis, Chao1, FB, and Shannon.
The dataframes within the two lists have same element names, but they contain different information. This is because list x has been generated after some operations on list my_list.
Now I want to write a function and run that on list x with `lapply. the command I want to write as a function is
forest(x$Chao1, slab = paste(my_list$Chao1$StudyID))

Here, the first Chao1 is coming from list x while the second  Chao1 is coming from the list my_list. It works for a single dataframe but I am unable to write the codes with lapply so that it runs across all the dataframes of the lists.
The desired outputs are four plots from four dataframes of list x.
Can anyone please help me?
Update: @gringer's approach works, but it saves all the plots in a single list. When I tried to save the images with following command it just saved the last image:
tiff(filename = "~/Desktop/image.tiff", width = 30, height = 30, units = "cm", res = 1200)

Is there a way to save all them individually?
Thanks

I am still quite confused. You mentioned that the names in my_list and x are the same, why do you need the names from my_list then? Also note you can [edit] your question and add all these details. Try for example, showing bits of the data.frames in the two lists, or something that would make us understand why you need two lists, not one.

Assuming the lists are named, and you want to fetch the element with the same name from each list, I'm guessing something like this would work (or at least get you closer to a solution):
v <- lapply(names(x), function(var){
   forest(x[[var]], slab=(my_list[[var]])$StudyID)
})

When creating images that don't have a multi-page format (i.e. anything except PDF), and not changing the image name with each save, you need to add in a format specifier to tell R what to do with new created images. For example, this will add a three digit zero-padded number for each new TIFF image:
tiff(filename = "~/Desktop/image_%03d.tiff", width = 30, height = 30, units = "cm", res = 1200)

Thanks gringer, it really worked. Being new in R, I would like to ask how does this line of code work? I mean, why I have to put double brackets [[]] there? Thanks


text:
How to compare Log2 Fold Change values
I have 5 sets of Log2 Fold Change values pulled off DualSeqDB. Each of the 5 groups corresponds to a different microbial taxon. Within the 5 groups, there are many Log2 Fold Change values, each corresponding to a gene that was significantly affected by the presence of the microbe.
My question is: is there any way to get the average Log2 Fold Change for each group? How could I calculate that?
The end goal is to compare each microbe's effect on gene expression, hopefully to prove that the phenotype of the microbe is in some way related to how much it can change the gene expression in its host.
I have considered using p-values to compare each taxa's effect on gene expression, but because all the datasets are significant, the p-values I got using fisher's method are all 0, so I don't know how I could go ahead and compare that to categorical variables (like presence of a certain phenotype).
Thanks so much to anyone who can help!


text:
Consensus sequence from vcf file
I am generating a consensus sequence for SARS-Cov-2. I am visually inspecting every candidate SNP in a pileup.
Sometimes I find that there are two consistent variants for the same position, but I can only generate one consensus sequence.
Which variant should I target:

The majority variant
The minoritarian one if I find that this SNP is under-reported in other public samples?
Are both required

What one is best practice?

Its a fundamental question. It would be useful to understand a little more about your sample

It was long ago... What do you need to know?

thanks, let me think a lot and I will formulate a question on the site and drop a comment to flag notify you

Making a consensus sequence tends to involve some arbitrary decisions. In this case, the common practices are:

Use the ancestral allele (this isn't always apparent)
Use the most common allele
Randomly choose

I think 1 and 2 are generally preferred, with 2 being easier to implement.


text:
Which sequencing technologies are considered short read
I couldn't find any information on this online. Which of the following sequencing technologies are considered short read technologies? Illumina, Oxford Nanopore, Ion Torrent, Pacific Biosciences

https://en.wikipedia.org/wiki/DNA_sequencing

Dear @bittu, this is the second question in a row that is nearly undoubtedly a homework question. It's alright to ask these, but you must disclose it!

As mentioned by Greg, list of short-read sequencing platforms is listed on wikipedia.
Of those you mentioned Illumina and Ion Torrent are short-read sequencing platforms.
Usually what makes short reads short is a PCR step in the library prep or directly on the flowcell to amplify the signal, which is the limitation for how long stretches of DNA are possible to sequence. Long read sequencing technologies sequence directly isolated single molecules.


text:
Which format is the most raw format
I wanted information regarding which format is considered the most raw from an illumina sequencer ?

fasta
fastq
bcl
bam
As per my research it should be bcl but I am not sure.

What do you mean by “raw”? None of the files you’ve listed are “raw” in any real sense. In fact, FASTQ and BCL are more or less equivalent, just stored in different ways. For the difference between FASTA, FASTQ and BAM, refer to https://bioinformatics.stackexchange.com/a/385/29.

BCL (binary base call) files are the most raw data that you will likely get. As a matter of fact, the images taken by the machine each cycle are even more raw but it makes no sense to interact with this data directly. The base calls are derived from counting the spots (clusters) on the images.
For any useful application you would convert them to fastq (e.g. using illumina bcl2fastq).
What is the goal of your question? Do you suspect some problems with your data?
Adding this information will make it easier to interpret what kind of answer you are looking for - and less likely that your question is downvoted.


text:
Batch effect correction using a subset of samples - using DE genes
I have some RNA-seq data with two very obvious batches as you can see in the PCA:

The samples of interest (A - H) are from tumor tissues. In addition, there are data for two cell lines I and J. The cell lines should be the same across the batches, while for the tumor samples there may be some additional biological difference.
My question is, if there is a way to use the cell line data to estimate the technical part of the batch information since there should very little biological difference here. Then I would like to use this to correct the technical aspect of the batch effect in the tumor samples.
When I run normal batch correction like ComBat or add the batch to the DE model the batch will be estimated over all samples including biological differences.
Update
I figured I could determine which genes are differentially expressed between the batches in the cell lines only. As mentioned before, the cell lines are highly similar biologically. Thus these DE genes should represent the technical aspect of the difference.
The plot below shows the PCA when after removing genes DE in the cell lines. Obviously, the differences for the cell lines are greatly reduced. Also the tumor samples are now less clearly separated.

Most of the removed genes were from shorter RNA species and a look at sequencing QC later confirmed that in fact the fragment size in the second batch was lower. Therefore I'm positive that I could identify and remove the technical component of the batch effect.
tl;dr
My data contains a strong difference between batches even for samples that are biologically the same (untreated cell lines). I looked at differentially expressed genes in the cell lines. This revealed to me that the second batch contained a lot more short RNAs due to lower fragment size of the library.
In the end I'm still wondering if there is an elegant solution for this requiring less footwork

The batch effect is expected to be the same across all samples - therefore you want to estimate its effect across all of the samples. Just add batch as an additional parameter to the linear model (this is presuming samples sequenced in the two batches have been randomised).
You can show the effect of "batch effect" removal with a PCA plot by plotting the residual error, e.g. in R if using DESeq:
mypca <-  prcomp(t(assay(varianceStabilizingTransformation(dds)))
pc.resid <- resid(aov(mypca$x~batch.info))
ggplot(pc.resid,aes(x=PC1,y=PC2) + geom_point()

(I haven't checked the code, but it will be something like this).

Thank you for the answer Greg. Unfortunately the samples are not randomly distributed across the batches. In the tumor samples there is a suspected biological effect on top of the technical batch (different regions of the tumor were sampled). Since this is not the case for the cell lines I was hoping to use them to estimate the "pure" technical component

Ah - unfortunate. Based on the PCA, I'm not certain I'd do anything other than treat the results with a bit more caution - your idea sounds reasonable though (if you can find a solid method)


text:
Combining read counts from three separate GEO studies
I want to do differential expression analysis with DESEQ2. I have three read counts files downloaded from GEO (small RNAseq based) where the number of miRNAs and id is nearly the same. These studies have used the same tool for alignment but different tools for counting aligned reads.
My question is

Can I combine these read count files
Is there any recommended meta-analysis based tool for this purpose
Is there any other way to do this other then redoing RNA-seq analysis.

Thank you

Yes, but there will be a batch effect.
Not that I know if, you'll need to use sva or something along those lines to handle the batch effect.
Crossing your fingers and hoping the batch effect is small is pretty much your only other option if the data isn't setup in a way that you can use sva or similar. So yeah, redoing the experiment may be the only real option.


text:
What is the file .mums created by mummer?
I am doing genome sequence alignment using MUMmer, in particular I want to do a dotplot with mummerplot. So the passages that I did are:
1.create a file .mums with the following command line:
mummer -mum -b -c H_pylori26695_Eslice.fasta H_pyloriJ99_Eslice.fasta > mummer.mums

2.and then plot all the MUMs with MUMmerplot:
mummerplot -x "[0,275287]" -y "[0,265111]" -postscript -p mummer mummer.mums

The example is taken here. The plot works, but I would like to understand several things to have more consciousness of what I am doing: the file .mums has three columns with many rows. What do they represent ? I know that it is a little explained in the link given before:

This command will find all maximal unique matches (-mum) between the reference and query on both the forward and reverse strands (-b) and report all the match positions relative to the forward strand (-c). Output is to stdout, so we will redirect it into a file named mummer.mums. This file lists all of the MUMs of the default length or greater between the two input sequences.

So the file should report all the Maximal Unique Matches between the two sequences but I do not understand what are these three columns...
Moreover, can someone explain in simple words the reason of the term "unique" in "Maximal Unique Matches" ?
Thank you in advance .

Maximal Unique Matches
To answer the question about MUMs, there are two important definitions:

A "match" is maximal if it cannot be extended in either direction and still be a match i.e, for two strings $s,r$, a match
$s_i,...,s_{i+k}$, $r_{j},...,r_{j+k}$ is maximal if $s_i,...,s_{i+k+1}$, $r_{j},...,r_{j+k+1}$ is not a match and $s_{i-1},...,s_{i+k}$, $r_{j-1},...,r_{j+k}$ is also not a match.
A "match" is unique if the matched substring is not present anywhere else in either genome. i.e, for two strings $s,r$, a match
$s_i,...,s_{i+k}$, $r_{j},...,r_{j+k}$ is unique if the substring $s_i,...,s_{i+k}$ is present only once in each genome.

With those definitions in hand, we can easily define a Maximal Unique Match (MUM) as an exact match that is both maximal and unique. Note that MUMmer also provides support for Maximal Exact Matches (MEMs), which do not require uniqueness of matches.
MUMmer output
As far as the output goes:

For each match, the three columns list the position in the reference sequence, the position in the query sequence, and the length of the match respectively. Reverse complemented query positions are reported relative to the reverse of the query sequence unless the -c option was used.

For example, the row:
1      1120       299

denotes a MUM of length 299 starting at the first base in the first input file and the 1120th base in the second input file. For matches on the reverse complement strand, you'll notice a second part of the output that contains the matches on the reverse complemented strand.
Here is a sample output I have on my device:
>$ $  mummer -mum -b -c England1.fna Indiana-USA-1_Saudi_Arabia_2014.fna
> gi|633896549|gb|KJ813439.1|
       1      1120       299
     301      1420       174
     476      1595        65
     542      1661        81
     624      1743        43
     668      1787       264
     933      2052      1042
    1976      3095      1299
    3276      4395       330
    3607      4726        69
    3677      4796       885
    4563      5682       812
    5376      6495       792
    6169      7288       831
    7001      8120       731
    7733      8852      1963
    9697     10816      1176
   10874     11993       617
   11492     12611        41
   11534     12653       164
   11708     12827      3007
   14716     15835      1040
   15757     16876       489
   16247     17366       412
   16660     17779       337
   16998     18117       463
   17462     18581       325
   17788     18907       686
   18475     19594       599
   19075     20194       342
   19418     20537       179
   19598     20717       532
   20131     21250       152
   20284     21403       935
   21220     22339       348
   21569     22688       764
   22334     23453       351
   22686     23805        73
   22760     23879        29
   22790     23909       180
   22971     24090       532
   23504     24623       143
   23648     24767       107
   23756     24875       434
   24191     25310        59
   24251     25370        47
   24299     25418       214
   24515     25634       224
   24740     25859      1386
   26127     27246        39
   26167     27286        63
   26231     27350       114
   26346     27465       102
   26449     27568      1399
   27849     28968       207
   28057     29176       938
   28996     30115       726
   29743     30868       369
> gi|633896549|gb|KJ813439.1| Reverse

Now, if you compute the reverse complement of the input and rerun, you'll notice the matches appear in reverse in the second section denoted by > gi|633896549|gb|KJ813439.1| Reverse

Thank you very much @Throckmorton . But I have still a doubt. How can mummerplot do the dotplot of the MUMs if in the .mums file contains only the "x starting value" and the "y ending value", let's say, (and the length) of the MUM ? Moreover, I do not understand why, using the options -b and -c (so running the first command line written in my question post), in my file there is also the list of the matches on the reverse complemented strand...

For each row $x,y,z$ in the mummer output, there will be a line between points $(x, y)$ and $(x+z, y+z)$ in the plot representing the MUM.

@Manuela There are matches on the reverse complement strand because you used the -b flag. If you run mummer -h you will see what the different flags do. For example, -b             compute forward and reverse complement matches

ok @Throckmorton for the first comment (I actually read your answer wrong, sorry). For the second I still can not understand because also you runned the same command with -b and -c and you do not have the list under "reverse"... isn'it it? And also in the manual (as I mentioned) there is: This command will find all maximal unique matches (-mum) between the reference and query on both the forward and reverse strands (-b) and report all the match positions relative to the forward strand (-c).

The reason I have no matches for the reverse complement is because the two sequences I ran MUMmer on are for the same strand and don't contain any inversions. If I invert the second half of the query sequence, those matches will then appear under the "Reverse" strand.

As far as the -c flag goes, consider the string ATTTATTTGC and the reverse complement GCAAATAAAT. If I have a match on the reverse complement strand at position 2, that's actually position 9 on the forward strand. With the -c flag, mummer will report the match at position 9 as opposed to reporting it at position 2 without the flag.

ok thank you very much !!!! @Throckmorton

found this in the manual, so I guess the first two columns are the reference and query positions and the third is the length. I am not 100% sure either, so let me know if you find out more.
http://mummer.sourceforge.net/manual/#mummerplot

Just a simple three column list, noting the position and length of every maximal exact match. Note that for reverse complement matches (produced with the -r option), the query start positions will reference the reverse complement of the query input sequence.

Thank you @HermannGregorDallinger ! Yes then I found an explanation like the yours always in the manual but in a different part of it that is:

For each match, the three columns list the position in the reference sequence, the position in the query sequence, and the length of the match respectively. Reverse complemented query positions are reported relative to the reverse of the query sequence unless the -c option was used. As was stated above the -L option adds the sequence lengths to the header line and the -s option adds the match strings to the output, if these options were used the format would be as follows: etc...

But I still wonder if these query and reference positions are the "starting" positions from the left or from the right... they should represent like a point (x,y) in a Cartesian plane but I do not understand very well which point.

Looks like this is undocumented, but a sensible default would be to get the end position by adding start and length, so they would be as you said from the "left" and the "right" position would be start + length. For most applications and for producing a "rough" dotplot, I guess it will not matter. You could try using blast to get more exact and detailed output if it is necessary.


text:
Automatize classical debugging in a bash script? (normal debugging by hand)
And thanks in advance. I am debugging a text file which is the one the two inputs on the program RGmatch https://bitbucket.org/pfurio/rgmatch/src/master/.
The code works. Is that some records give division by zero.
I am doing it by hand.
But is there a way by making a bash script, to see which records show the problem?
I do not want to touch the python code.
python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b ../../test.bed -o myassoSort.txt

I have tried with classical debugging, but my actual question is how to automatize that?
Also there are no empty records. The bed file was created based on refseq.txt from the latest human assembly. And also the GTF file. On the refseq.txt I processed to obtain a bed file.
I have checked if this works on a single record, it does, but not in some specific ones. The error I get is:
  File "rgmatch.py", line 1311, in <module>
    main()
  File "rgmatch.py", line 287, in main
    run(gtf, dhs, outputfile)
  File "rgmatch.py", line 1153, in run
    pctg_area      = (float(region_overlap)/intron_length)*100
ZeroDivisionError: float division by zero

The record I found so far is this one.
chr19    14087850    14090751
I did grep on both GTF files and Refseq Files and in both files there is information about the record. Because the program search for areas between region of interest, If I put inside just the record only in GTF file and bed Files no error occurs.
grep "14090751" refseq.bed
692  NM_138352.3 chr19   -   14087850    14090751    14088416    14090420    6   14087850,14088670,14088879,14089067,14089784,14090084,  14088575,14088782,14088976,14089362,14090084,14090751,  0   SAMD1   cmpl    cmpl    0,2,1,0,0,0,
grep "14090751" hg38.ncbiRefSeq.gtf
chr19   ncbiRefSeq.2021-02-10   exon    14090085    14090751    .   -   .   gene_id "SAMD1"; transcript_id "NM_138352.3"; exon_number "1"; exon_id "NM_138352.3.1"; gene_name "SAMD1";
chr19   ncbiRefSeq.2021-02-10   5UTR    14090421    14090751    .   -   .   gene_id "SAMD1"; transcript_id "NM_138352.3"; exon_number "1"; exon_id "NM_138352.3.1"; gene_name "SAMD1";

Please include the things I asked for over at [linux.se]. We need to see the actual error message. Also, I call it "classical debugging" but I don't think anyone else will understand the term :) So, please [edit] your question, add the specific error message and then show us the commands that you want to automate. Have you checked it it works with just a single record? Does it work with just the first line of the bed file?

Also, if you still need help, feel free to ping me in the site chat room where we can directly debug the issue.

If you want to just start incrementally adding lines from the bed file until you hit the one that reproduces the problem, you can do something like this:
for i in $(seq 20 30 $(wc -l < ../../test.bed)); do
    head -n $i ../../test.bed > slice.bed
    if python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b slice.bed -o slice.out
    then
        true
    else
        echo "Failed when running with the $i first lines."
        break
    fi
done

Explanation
seq 20 30 $(wc -l < ../../test.bed) will iterate over all numbers between 20 and the number of lines in your bed file with increments of 30. Next, the head -n $i ../../test.bed > slice.bed will take the first $i lines of the bed file. This way, the loop will run the rgmatch command on the first 20 lines of the file, then the first 50, then the first 80 and so on as long as the command runs successfully. When the rgmatch command fails, it will report the value of $i to help you find the line that caused the problem.

Thank you @terdon I will edit the answer and try this.. If not I will appreciate very much your help. Right now I have a work meeting. But after that, I will do everything you said to me. Thank you !

De res, @EdmondGeraudAguilar! And do feel free to ping me in chat if you need more help.

It works ! I mean, I would that the script makes like a try-catch then store the results.. but I will handle with this. Than you very much

@EdmondGeraudAguilar it's very easy to extend to whatever you need. But this is why I asked you to edit and explain what commands you want to automate. If you want to keep the results of each run, just change the rgmatch output file: if python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b slice.bed -o slice.$i.out; then ...

For those who may use this function be aware that some start coordinates of transcripts are missing in the GTF files that is the problem


text:
How to specify the list of contigs in AGP gap file used for ENA submissions
I am battling the submission of scaffolded assembly to EBI, somehow I can not get AGP gap file right.
background
The only option seems to be using the Webin-CLI submission interface. For the scaffolded assembly submission, one needs the manifest file, the sequence (fasta or flat file), and AGP file, which denotes how contigs are glued together into scaffolds. I got already that the sequence file should contain only contigs, and EBI will glue them together using the AGP file.
the problem
So far so good, I checked the specification and the example provided and generated similarly looking files for my assembly, but then when I tried to validate the files using
ena-webin-cli -context genome -userName XXX -password XXX -manifest manifest_file -outputDir asm01_submission -inputDir ena_submission_input -ascp -validate

I got tons of errors with type
ERROR: Scaffold "NODE_667656" has only "1" component, minimum two components expected.
ERROR: Scaffold "NODE_353389" has only "1" component, minimum two components expected.
...

supposedly for each contig that actually was not scaffolded. I am quite confused about this, as the example file contains tons of scaffold with one component only (for example the first one EG1_scaffold1).
What am I missing?
Sample of my AGP file (incl 1 scaffold and 3 contigs):
Path_2003       1       1672    1       W       NODE_32694_length_1672_cov_36.024451    1       1672    +
Path_2003       1673    1772    2       U       100     scaffold        yes     align_trnscpt
Path_2003       1773    2521    3       W       NODE_61726_length_749_cov_25.233631     1       749     +
Path_2003       2522    2621    4       U       100     scaffold        yes     align_trnscpt
Path_2003       2622    5705    5       W       NODE_20424_length_3084_cov_26.290988    1       3084    +
NODE_667656     1       86      1       W       NODE_667656_length_86_cov_136.111111    1       86      +
NODE_353389     1       155     1       W       NODE_353389_length_155_cov_25.230769    1       155     +

Oh you silly sausage! Don't use AGP gap files.
This is not a direct answer but a workaround. However, if your scaffolded genome is annotated, you might as well submit a flat file (.embl) instead of fasta and the gap file. The thing is, making that file is actually a lot more straightforward thanks to amazing EMBLmyGFF3. You will also need the locus tag (a string used as prefix for loci in your genome) you must pre-register in the web interface (it's among "Project" options). Note that although the script does handle some EBI specifics (like automatic filtering of gff3 tags they don't like, or filtering sequences < 100 bases), some there are things that need to be taken care of. For example, I had to remove duplicates (using agat_sp_fix_features_locations_duplicated.pl) and overlapping genes (using agat_sp_fix_overlaping_genes.pl), both available in AGAT collection. This might be different for genomes annotated with a different tool (mine was braker2).
Anyway, once I was advised to use this tool, it took like 2 hours for 400Mbp genome including the uploading time, so I certainly recommend it over wrestling with AGP files.


text:
How to interpret complete genome sequence in NCBI
I have a very basic and quick question: does the complete sequence of nucleotides of an organism (for example this one that I am considering in my homework of Pyrococcus Fuoriosus) reported in NCBI represent the sequence of nucleotides of the forward strand (5' to 3') of DNA?


text:
How to graph sequence ranges, for example intron sequences or conserved domains
I am a bioinformatics beginner.
I want to draw line graphs such that for a gene, say gene_A, intronic regions are (20,100), (500,600) and a conserved domain was found in (400-550), then I should get base position/sequence length on the x axis and gene name on the y axis.
Is there any software that can help me do this?
I used matplotlib but would like better visuals, any suggestions would be helpful. Thanks!

Please [edit] your question and give us some more details. What will your input be? What exactly are you trying to graph? You say "gene name on the y axis" but show an example with fractional numbers instead. You say "base position/sequence length on the x axis" but those two are very different things. You say you have an intron from position 500 to 600 but then have a conserved domain at 400-550, so part of the intron. None of this makes much sense, I'm afraid.

Would ggbio do the trick?


text:
Label and labelbox mismatched in Pymol after ray
I would like to ask how to overcome the problem of mismatched boxes and labels (pic). Is there an option focused on this?
thanks
Radim

It is a glitch in PyMOL.
Caveat. I would say that, whereas cmd.label is great for adding labels en-mass for internal figures, it is not great for figures for dissemination (which require few strategically placed labels possibly with a faint white outer glow —cf. your D330): most figures in papers I would say are Powerpoint or Photoshop labelled.
You called it 'labelbox', this is not a normal term, so you may be using a custom script. The box is normally called a bounding box and connector line. To add it you need to set a few values. You may have done the following:
label name CA and resi 44, 'Foo'
set label_outline_color, 0xff8c69  # font colour to salmon
set label_connector_color, 0xff8c69 # connect or bounding box
set label_bg_outline, true  # adds a border (lowercase bool)
set label_connector, true   # adds a line
set label_connector_mode, 1 # how the line looks label_position
set label_position, [5, 5, 5]  # offset the label
set ray_label_connector_flat, true # cylinders as line? No thanks.

The issue is the set label_position, [5, 5, 5]. If you give the camera-based reference vector [0,0,5] you get no line as Z is depth. This ray traces fine. Otherwise the glitch you found happens. My guess is that the box (but not the label) gets calculated based on the viewport window size, not the ray trace window size. So my guess is that there is a resolution and ray trace window size combination just right that they do align... but on my machine at least this is low res.
To make matters worse, the shadow is horrendous and in a random location and label_shadow_mode does not stop it. Let's not even discuss ray_label_connector_flat. So the best way to avoid this is to simply label in Photoshop, which is much more convenient anyway.


text:
How to compute a p-value for the similarity measurement of two DNA sequences/contigs?
I'm attempting to come up with a statistical model to determine the p-value for the similarity value $ 0 \le s \le 1 $ of two DNA contigs.
The idea is to get a numerical measure to estimate if two contigs are homologous.
I want to account for base composition (and/or other structural features) if possible, but not necessarily because I might want to do the analysis only from the information about the length of the sequences/contigs and largest overlaps. Maybe something inspired from other statistical significance in bioinformatics .
Assume we have two DNA sequences: $a$ and $b$ with the basepair lengths $a_{\text{len}}$, $b_{\text{len}}$. Their overlap $o_{a-b}$ is the number of contiguous bases that appears in both sequences. If we define the similarity $s = \frac{o_{a-b}}{\min{(a_{\text{len}}, b_{\text{len}})}}$ how can we define a p-value for $s$?

Okay, can you please give some more context in your question about a specific use case for this algorithm? It will help to get better answers.

What is basis of your denominator ? How will you calculate ? what does  ,  comma here, stand for? How is it better measure ?

Well, here my thoughts:
The equation $s$ somehow tries to represent the ratio $\frac{equality}{equality+inequality}$. But $s$ returns only a value. And to assing the p-value to a single value... is complicated.
How could you approach this?

Using some kind of Hamming distance for each base pair at the position $n$ ($1 \leqslant n \leqslant len$) between the sequence i ($s_{i}$) and sequence j ($s_{j}$) in order to get a set of binary data:

$$s_{ij} = \begin{Bmatrix} 1\ if\ seq_{i_{n}}=seq_{j_{n}} \\ 0\ if\ seq_{i_{n}}\neq seq_{j_{n}}   \end{Bmatrix}$$
With this, you get a set of binomial data (eg $\{0,1,0,1,0,1,1,1,0,1,1,1,0,...\}$). Now, you should be able to perform some (Binary?) test and get the p-value.
PD: Consider aligning the two sequences if there is interest in discussing base pair homology.

The overlap is what the aligner labeled as the largest similar subsequence

Maybe something along the lines of how blast compute p-value


text:
Conversion of SAM to BAM files
I am very new to micro RNA analysis. I have been using H. sapiens, GRCh38 + major index as given in the Bowtie Website to align with my trimmed FASTQ file .
The command I am using to make very sensitive alignments with Bowtie2 is given below.
module load  Bowtie2;
bowtie2 --local -p 8 -q --phred33 -D 20 -R 3 -N 0 -L 8 \
    -i S,1,0.50 -x grch38_1kgmaj \
    -U /gpfs/ycga/scratch60/nicoli/ag2646/HumanmicroRNAbowtie_verysensitive_1.06.2021./3kPa_HDF_miRNA_1bbduk.fastq \
    -S /gpfs/ycga/scratch60/nicoli/ag2646/HumanmicroRNAbowtie_verysensitive_1.06.2021./3kPa_HDF_miRNA_1bbduk.sam

Then I am using:
samtools sort -bS 3kPa_HDF_miRNA_1bbduk.sam > 3kPa_HDF_miRNA_1bbduk.bam

to convert sam files into bam files.
However, I am hitting the error described below:
[W::sam_read1_sam] Parse error at line 46371

When I grepped the specific line it shows:
SRR8248790.46371        4       *       0       0       *       *       0       0       CCTGGATGATGATAAGCAAATGCTGACTGAACATGAAGGTCTTAATTAGCTCTAACTGACTAGATCGGAAGAGCAC   BBADAGEHIFHHIIIIIGGFIHIIIIIIIHHHHHIHHHEHEHEHFEHIIIEHIIIIIEHIIHIIIHGHHHIIII1C   YT:Z:UU

Any specific help troubleshooting this error will be very useful.

The [W: is just a warning. If the process errored out, there should be another actual error message that follows that warning.

Is it the last line of the samfile? I wonder if the sam file could be truncated (like a missing newline or something).

How did you "grep the line"? Is that actually line number 46371 of the file, or is that a random line that happens to contain the string 46371?

Is your samtools sort command actually a samtools view command? The -bS flag would not work with samtools sort: http://www.htslib.org/doc/samtools-sort.html. I tried it out of curiosity and got samtools: invalid option -- b.

See similar situation here: https://www.biostars.org/p/169277/. Have you confirmed that bowtie2 finished without errors? Does your sam file look ok? how many lines does it have?

My bowtie file alignment completed without truncation but only showed errros.

Are you able to trace the same line in the raw fastq file (before trimming) with the same identifier?? If yes, what is the content of that line?

I have resolved all the errors. Thanks for all your help.