text:
How important are the homozygous variants that get unnecessarily deleted using liftover?
I'm referring to the text described here:

These tools [NCBI remap, CrossMap] operate only on the sites present in an input VCF, and return the representation of those sites in a new genome assembly. This does not capture all variation, however. Consider an individual who has sequence reads that indicate they match the GRCh37 reference genome assembly at position GRCh37.chr1:169,519,049 (i.e. the individual's genotype is T/T). Because the individual is homozygous reference at that site, there will be no variation present in their VCF file created on GRCh37. However, the analogous position on the updated GRCh38 reference genome assembly, position GRCh38.chr1:169,549,811, has the reference base C. Consequently, if the individual's read data were analyzed on GRCh38, they would be identified as homozygous for a C->T SNP. Because this site is not present in the input GRCh37 VCF, it is never added when creating a GRCh38 VCF by these other tools.

How important is this? Do you think it is relevant to any types of analysis?

I'm mostly just wondering if the bioinformatics community thinks the added benefit is large enough to be relevant. I haven't encountered a need for this in a project.

All this is saying is that the tool doesn't lift over variants which aren't present in the vcf. I can't imagine this being a big deal in any analysis I've done and it's exactly what I'd expect the tool to do. There's no rule that invariant sites (with respect to the reference) can't be in a vcf, so if you are concerned then just retain the variants in the vcf file at the calling stage


text:
How do I figure out pangenome location of hg38 coordinates?
I want to a mapping that has linear coordinates and corresponding pangenome (e.g., .gfa format) locations. For example, this type of information could say that region chr 1 1000-2000 corresponds to nodes 104-108.
I have tried vg find, vg indexing tools, and manual .gfa analysis, but I could not get this information.
Any of these would be a solution:

A description of an algorithm to extract this information from a .gfa file
A file where this mapping is explicitly given
A tool that works to retrieve this information

The keyword to look for is 'liftover'. Here's liftover information for translating from GRCh38 to CHM13:
https://github.com/pangenome/chm13-grch38-liftover
If this isn't what you need, I suggest asking the pangenome people directly via a github issue:
https://github.com/pangenome/chm13-grch38-liftover/issues

source: Stack Exchange

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

sort: votes


text:
Is there a database or file containing the type of variant (e.g., if (non-)synonymous, in exon or not) for all human genomic positions?
I'm looking for data which shows the type of SNV variant (at least if the variant synonymous or non-synonymous, and ideally if the variant is in a CDS or not as well) for every variant with respect to the hg38 reference.
For example, a table that shows something along these lines for all variants:
chr1, position 1234, genotype T, reference genotype, not in CDS

chr1, position 1234, genotype A, synonymous variant, not in CDS

chr5, position 1111, genotype G, non-synonymous variant, in CDS

I don't care which format the data is in or how it is laid out as long as some of this information can be extracted. I don't have a preference between a database of all possible variants (including variation which is not observed in reality) and a database of all variants (polymorphic sites), from, e.g., 1000 Genomes or a similar panel.
I am not looking for deletions, insertions, or CNVs (though it is fine if the dataset has them; I will just exclude them). I'm looking for the annotation (at least if it is synonymous vs. non-synonymous) for every single base-pair variant.

When you say for every variant, I assume you mean for "every single-base missense variant", right? No deletions, no insertions, no multi-nt variants. Please clarify that.

@terdon Yes, I'm looking for all single nucleotide variants. It doesn't have to include SVs, CNVs, etc.

I think the CADD database has this information https://cadd.gs.washington.edu/download. You want 'All possible SNVs of GRCh38/hg38 incl. all annotations'.

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


text:
Where do I get a large reference VCF?
I would like to download a large .vcf file containing many (hundreds or thousands) of samples. Ideally, I would download different population-specific .vcf files, but the ability to sort/filter by ancestry group is fine. Where do I get such a file? I prefer GRCh37 for consistency with other files I'm using.

Update: I tried the HGDP file--I'm trying the 100 genomes one now, in addition to trying http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/release/20190312_biallelic_SNV_and_INDEL/ now.
I downloaded the file, but when I opened it it seemed to have no samples (instead of having many samples). Here is a full row:
 1  10145   .   AAC A   151.43  PASS    AC=4;AF=0.0127389;AN=314;BaseQRankSum=0.72;ClippingRankSum=0;DP=24144;ExcessHet=160.275;FS=0;InbreedingCoeff=-0.1558;MLEAC=5;MLEAF=0.002604;MQ=1.63;MQRankSum=0.764;NEGATIVE_TRAIN_SITE;QD=3.52;ReadPosRankSum=0.085;SOR=0.458;VQSLOD=-0.9391;culprit=QD;VQSRMODE=INDEL;NS=157;ExcHet=0

This question has also been asked on Biostars

@BigMistake please [edit] your question and give us some context on what you need. You have now mentioned you also need "ancestry information" (which isn't very clear, what kind of information?), and of course you need samples, but please explain all this in your question so we can know what you are after. Ideally, tell us something about how you will use the file so we will be able to give you better answers.

Hi @BigMistake thank you for posting your answer. Could you kindly mark this as "accepted"?

Solution: Use https://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase3_v5a/b37.vcf/, then to get the ancestry information (i.e., the ancestry group of each sample, e.g., "HG00566" or "NA19758") use https://www.internationalgenome.org/data-portal/sample.

For anyone who arrives here looking for a large reference VCF but not wanting GRCh37, you might consider the HGDP+1KG subset of the Genome Aggregation Database's dataset. You can find a VCF file per chromosome aligned, I believe, to GRCh38, here.


text:
Does pooled sample genotyping microarray analysis exist?
I want to genotype samples in small batches with a common microarray (e.g., the GSA). However, the cost is very high for small batches.
With a higher sample throughput/volume, per sample cost becomes more reasonable as the processing cost for a single run is spread out over more samples.
I'd like to process no more than 16 samples per run, ideally even one or four samples per run. (Being able to submit a single sample for microarray analysis without having to pay to run the rest of the empty chip would be great.)
It is possible that some company or core offers "pooled" microarray services, where I could submit a small number of samples that would be processed with other groups' samples on the same microarray run. This would mean there are more samples per run, so it could potentially cost much less.
Looking on core and company websites, I don't see this offered anywhere.
Does this service exist?

I understand the question entirely. My personal view is this is primarily a wet-lab question. Companies will do mass chip printing to reduce costs, which are stored. Grouping between independent clients I would suggest is avoided because of the risk of any form of "spill-over". For example, chips get mislabelled and you get someone else's results.  For the company it's a disaster and they lose client confidence. Better to work alongside colleagues with similar aims as you for the "pool".

The experimental problem with microarrays/chip technology is the variation in quality. Statistically this can be corrected, but wet-lab perspective is they don't like any inclination of sub-optimal quality.

@M__ Yes, that makes sense. I imagine that's why most cores don't operate like that

This can also be done at the sender (rather than the chip runner) level. Multiple senders get together to send samples


text:
Salmon Pseudo count when dealing with male and female RNA-seq data
I've generated a quant seq data that I intend to use to compare male and female gene expresion, with a focus on sexual chromosome.
For my species (three-spined stickleback), it is a classic XY sex determination system, and i have a male and a female reference transcriptome, including or not the Y chromosome.
I am not sure how to handle that with SALMON, and I have 2 questions in particular:

Given that I use two different transcriptomes depending of the sex, can I still compare the obtained normalized pseudocount of would it affect the estimation too much?
I want to compare expression levels of genes present in both the X and Y reference (i.e genes that are not recombining but have not been lost or degenerated on the Y). For that, I'd need to add the gene counts of the Y and X copy. Is it Okay to add  pseudo count from different gene copies, or will that be issue? I do not understand the pseudo-count statistical model enough to have an idea about that

My other option is to go with a STAR / HTseq pipeline, but my mapping rates seems about 20% lower with this tools.


text:
Mapping statistics from the bam files
I would like to find the mapping statistics from the sorted bam files.
Samtools flagstat gives the output only for a single file. What is the easiest way to find the mapping statistics for all sequences? How do I find the mapping percentage for all sequences? Is it possible to use the samtools flagstat if I merge all the files in to a single bamfile?

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

Hi @Kam and welcome. What sort of stats are you thinking of? E.g. coverage, average read length? Also what sort of sequence data is this?

I would like to find the mapping percentage for all the reads after mapping. I have around 300 bam files after mapping. Is it possible to find the mapping percentage for all the reads?

run it in a loop. something like for f in bamfiles; do ...; done

Which mapper are you using? Have you checked the logs of your mapper to see if these metrics are already provided there? Are you sure you are retaining unmapped reads in the final BAMs in order to calculate such rate?


text:
Create bed file where each column is the number of intersections with another bed file
For example I have 3 bed files:
window.bed:
chr1    0   1000
chr1    1000    2000

lp1.bed:
chr1    20  50
chr1    30  40
chr1    1500    1900

lp2.bed:
chr1    68  786
chr1    899 987

How can I create this output:
chr1    0   1000    2   2
chr1    1000    2000    1   0

In this output columns 1, 2, 3 are the standard bed format cols, of the chromosome, start, and end of the features taken from window.bed above. Column 4 is the number of intersections between window.bed and lp1.bed per feature in window.bed. For example, the first element in column 4 is '2' because there are two intersections between window.bed and lp1.bed. Column 5 is the same but for the number of intersections between window.bed and lp2.bed for each feature in window.bed.
In the case that a feature is a subset of another feature in either lp1.bed, or lp2.bed, they should be counted as two intersections. For example, the second feature of lp1.bed is a subset of the first feature. This is reflected in the first element of column 4 of the output.
I suspect that it is doable with bedtools intersect, but I cannot work out the syntax.

Via bedmap --count, which gives the option to specify a custom value when regions do not map:
paste <(bedmap --echo --count --unmapped-val 0 --delim '\t' window.bed lp1.bed) <(bedmap --count --unmapped-val 0 window.bed lp2.bed) > answer.bed

Ref. https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#element-and-overlap-statistics

The relevant bedtools option seems to be -c:
    -c  For each entry in A, report the number of overlaps with B.
        - Reports 0 for A entries that have no overlap with B.
        - Overlaps restricted by -f, -F, -r, and -s.

We can get the number of overlaps per region per file with:
$ bedtools intersect -c -a window.bed -b lp1.bed
chr1    0   1000    2
chr1    1000    2000    1

$ bedtools intersect -c -a window.bed -b lp2.bed
chr1    0   1000    2
chr1    1000    2000    0

So we're almost there, we just need to combine them, and that's a job for paste:

$ paste  <(bedtools intersect -c -a window.bed -b lp2.bed) \
         <(bedtools intersect -c -a window.bed -b lp1.bed)
chr1    0   1000    2   chr1    0   1000    2
chr1    1000    2000    0   chr1    1000    2000    1

But that gives us extra fields, so we can use cut to remove them:
$ paste  <(bedtools intersect -c -a window.bed -b lp1.bed) \
         <(bedtools intersect -c -a window.bed -b lp2.bed | cut -f 4)
chr1    0   1000    2   2
chr1    1000    2000    1   0

The paste command just joins lines from its two input files in the order it receives them, but the bedtools commands will always output the same order as the lines in window.bed so this should be fine.
The <(command) construct is called "process substitution" and is a way of passing the output of a command as a file to another command that expects files as input. Note that it isn't supported by all shells, but if you're using bash or zsh which are the most common these days, it should work.


text:
Is there any way to handle CONECT records in AutoDock Vina?
I modelled a macromolecule/receptor using Avogadro and exported it as receptor.pdb. Many entries in this file are CONECT records. If I try to run AutoDock Vina with it, this happens:

So, if I prepare this receptor again for docking, but, this time, removing the CONECT records (as I've seen many people do), my receptor file now exhibits broken bonds. Like this:

I've tried converting the .pdb file to .pdbqt with OpenBabel — removing all but the largest contiguous fragment and ignoring CONECT records as well —, but nothing works. I've also looked into it, and it looks like Vina doesn't have any parameters to handle CONECT records. Isn't there any way at all to do this?
Thank you in advance.

Hey, @M__! I'm sorry for not doing this sooner. I was a little short on time yesterday. What about now? Is there anything else I should add?

Hi, thats its better. If you copy and paste the ascii characters thats preferred to screen shots for the error output.


text:
2D diagram of intermolecular interactions in R
I've docked a large inorganic compound (receptor) with an organic molecule that isn't a protein (ligand). I wanted to see the interaction between them in a 2D diagram, but BIOVIA Discovery Studio (DS) doesn't work well with pdb files converted from pdbqt. And since the receptor isn't a protein or a fragment of DNA/RNA, DS also won't show the diagrams (the same with PLIP and Proteins Plus). I've also tried LigPlot, LeView and PoseView, but none of them worked as well. Does anyone know another  software that draws 2D diagrams, and, if not, is there any way to do this in R???

Hi @ginn are you sure R is the right platform for this job? Generally, this is where Python excels.

Oh... I hadn't thought of that yet  but thank you so much for your input. Would you know, perhaps, of any package that does this job...?

Beyond Bio.pdb https://biopython.org/docs/1.75/api/Bio.PDB.html and of course pymol - I dunno - I'm not a specialist in this area. There are members here who are experts, so maybe they'll respond.

I've tried PyMOL, but apparently you have to install the pose view plugin, which you need to be licensed to use...  But thank you so much once again!!! I'll try bio.pdb in the meantime ☺


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

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

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


text:
Reading same BAM file twice with htslib
I would like to read iteratively the same BAM file twice using htslib.
My option is to use hts_open, sam_hdr_read, then sam_read1 as much as needed to read the file until EOF.
Is there a possibility to "rewind", and use sam_read1 again from the beginning of the file, without closing the file, destroying the object, and opening it again?

EDIT
My aim is group the reads by the barcode.  The first pass would allocate the specific sizes first (the number of reads by barcode), then store the reads.
I know I could sort the reads first (using samtools sort -t BX), but I would like to skip this step using pre-allocation.

@M__ Thanks, I edited the question.  I hope it helps!

The best solution is to reopen the file. For a coordinate sorted file you have the option to jump back to the start of the first chromosome with index but that is more complex, less flexible and won't be faster.


text:
How do I calculate differential expression for RNA-seq values with the "limma" package and the "ebayes" function?
So for context, I have a set of TPM values for multiple genes for different samples, and I need to calculate the differential expression for RNA-seq values. I've been using this website as a guide: https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2019_March_UCSF_mRNAseq_Workshop/master/differential_expression/DE_Analysis.Rmd. The overall aim of this calculation is so that I can use PANDA to map TM motifs/combine gene expression data to create networks representing interactions between transcription factors and genes, and GSEA to analyze genes ranked by fold change or differential expression p-value.
So what I need to do is match one male to one female of the same age group (which I already did prior; the number of males outnumbered the number of females so I filtered the number of males so that there would be an equal number, and of the same age distribution). But it seems that I need to derive two factors and create a new variable "group" that combines factor1 and factor2. In this case, it would be the gender I would use as the group factors, right?

My advice is open a separate question on the topic of calculating DEG cut-offs on this site. I personally cannot  answer the first question (I don't know the project, specifically what the DEG is measuring and how pairwise DEGs are avoided), but I can present the rationale of calculating DEG cut-offs (i.e. different approaches used and why). Again this is the second part of the question and the SE (Stackexchange) approach is one technical request per question. In this instance the answer  is bit complicated (if stats aren't you thing) and is not a one line response in comments.

Okay, but do you have a recommendation on the specific number I use as a cutoff?

Which is..? It wasn’t clear to me if you had a number in mind or you were telling me to just eyeball it.

And how do you derive this?

Okay, I did what you asked


text:
Mapping single cell data with annotation
Entry level Bioinformatician here.
I am trying to map my single cell data with standard annotation file or assign my Seurat cluster to cell types. I have used the singleR and celldex package to map it but my problem is the single cell only consist of endothelial cells. I haven't found paper that only consist of rats endothelial cells. The single cell data is for rats.  I am stuck. What is the best thing that could be done here?

Given that it sounds like you already have a Seurat object, can you please show a UMAP of the cells that is coloured based on clusters identified in FindClusters? This will help us to understand the heterogeneity of the cells you are looking at.

It seems like you can annotate all your cells as endothelial cells, problem solved ;) But seriously, not clear what your expected cell types are (subtypes of endothelial cells perhaps) or the tissue type. If you can't find a suitable reference dataset for the reference mapping approach, then you can always resort to the 'old fashioned' way of unsupervised clustering (e.g leiden) plus checking marker gene expression to annotate the clusters.
An alternative approach, could be finding a suitable human/mouse reference and  then converting the gene ids to rat (ideally via 1:1 ortholog mapping) but this could get messy, so I suggest the former simpler way, at least as a starting point

Thanks, this seems to clarify the question. Is this what you mean @appleball?

Hello. thanks for the comments and concerns.
For sequencing the data, only endothelial cells were selected so the expected cell types would be the subtypes of endothelial cells, venous capillary or arterial. I have not had any luck in finding the annotation file for the endothelial cells.

Note that many cell selection methods are imperfect, so for your first pass you should start with confirming their identity by looking at the expression of markers of endothelial cells, and markers of other cells types (pericytes, muscle, ...).


text:
Phylogenetic tree from a large FASTA file
I want to create a phylogenetic tree using a 3.8MB FASTA file with ambiguous code. How can I perform multiple sequence alignment and the construction of phylogenetic tree with these ambiguous codes?
Could you please suggest some programs or R packages to complete this task?
How does MEGA11 software handle ambiguous codes? Will the alignment consider them to be mismatches?

Welcome @user17615. Firstly, I would recommend starting with a pilot study and then scale up. What is the genetic distances or taxonomic groups between the sequences, e.g. species, genera, families, orders (thats quite important)? Is it essential to do this in R, or what the requirements for R, e.g. is that the language you work in? Finally, how is your bash/Linux command line?

MegaX is poor for this purpose, it handles sequence ambiguities. The file sizes for complete phylogenetic analysis within R is going to be very tough, but its doable using "soft statistical approaches" for phylogenetic methods and can be processed via the ape library.