cannin/biostars_qa huggingface dataset:
Question: I intend to find common peaks ( in ChIP-Seq data, ) defined as peaks that exist in any pair of the 3 replicates I have. Of course there are 3 pairs, so peaks that overlap in 1 base in any pair of those would be considered valid.

If peak A and peak B overlap between two replicates, I will define the final "peak" as the union between the area that peak A spans and the area that peak B spans.

The issue is that I would like to have some way to **rank** the final peaks, at least to decide which peaks are "the best peaks".

Each final peak is a combination of *at least* two peaks; while the original peaks have meaningful p-values and q-values, I don't think their combination would be meaningful for the combined peak.

Can you suggest a meaningful way to rank the combined peaks?

Answer: You would need a statistical framework. I could think of several options:

1) Use a dedicated replicate-aware peak caller:

- [PePr][1], easy to use, accepts replicates, returns statistics per peaks which (from what I understand) respect the variability between replicates. It needs an input sample though.
- [Genrich][2], actually developed for paired-end data but can tweaked to accept single-end reads if you know the average insert size. Needs no input explicitely, combines statistics per peak using Fisher's method. Works in my hands well for ATAC-seq, did not test it for ChIP-seq. Method is currently unpublished.
- [CHIP-R][3], very new method I stumbled over on [Twitter][4], seems to be a pipeline that somehow assesses replicates in terms of reproducibility, had no time yet to look into it though, feel free to post feedback if you try it.

2) Naive intersections:

- Intersect a) with b) and the output with c), take the overlaps. That is one the one hand stringent as peaks have to be present in all three sets but relaxed as it does not consider the variability, so only checks for binary presence. Once you have this call peaks on the merged bam files, intersect with this peak list you just created, and then rank by the p-values (or any other metric from the narrowPeak file) for those called peaks that intersect with the intersection set. Probably the weakest of all the methods I describe here, it is rather a thinking-aloud, probably not very reliable, and quite tedious.

3) Use IDR:

- The Irreproducible Discovery Framework (here a [tutorial][5]) was developed to check for consistency in ChIP-seq data, but only accepts n=2. So you could take the two replicates with the best data quality and then use them to create a confidence peak set.

4) A rank-based meta-analysis.

- https://cran.r-project.org/web/packages/RobustRankAggreg/index.html might be an option. Rank the peaks of each sample by a metric of choice, then use this package to prioritize them based on rank-consistency.

As always it depends: What is the final analysis goal?

[1]: https://github.com/shawnzhangyx/PePr
[2]: https://github.com/jsh58/Genrich
[3]: https://www.biorxiv.org/content/10.1101/2020.11.24.396960v1
[4]: https://twitter.com/rhysnewell/status/1332127898854060033
[5]: https://hbctraining.github.io/Intro-to-ChIPseq/lessons/07_handling-replicates-idr.html
Question: I have whole genome sequences of several clinical isolates and I want to get the nucleotide frequencies (% A,C,T,G) at a given position. Among novel SNPs we could discover, there are some very specific positions that I want to look at. There are enough clinical isolates that I don't wish to this manually using IGV.

I can do this, but it's extremely inelegant and I think that there is some more elegant way that I am missing. I have summarize that there are probably two approaches.

**Approach 1**

1 . Run the basic samtools SNP pipeline, outputting all positions without a VCF ( remove -V) and outputting all possible reference alleles.
2. Write a script that parses the info column in the VCF file, get the depth info for each nucleotide, calculate the frequencies (for my variants of interest only); luckily I already have this script.

**Approach 2**

1. Use bedtools and convert that bam to fastq
2. Use seqtk and convert fastq to fasta
3. Use bedtools nuc command, specifying position of interest

I feel like this is too common a task to require the use of so many tools and file format conversions and even custom scripts, so I ask : what's a better way to do this?

Answer: What about parsing the relevant line from samtools mpileup?
Question: Hi,

I have 40 whole genome sequencing .bam files each one about to 1 TB. Lab manager has downloaded them from Cambridge university to her hard derive. In the other hand my computer in office is windows OS so I have to use linux OS in our Compute Cluster in university to analyze my data. For that I have to drag and drop each .bam file to my scratch in Compute Cluster that each file takes about 2 days for transferring and after finishing I am seeing .bam file becomes crashed. Now, it is about a month I have been given this data but I am still struggling with transferring them. If one of you were in my place what would you do?

I thought to ask lab manager to download these .bam files directly to Compute Cluster (but she has not done yet)
I thought to ask IT service to install linux on my computer although I guess I would need cluster computing again (

I really don't know what to do

Any consultant please?

Answer: That sounds like some big amount of data and therefore it will take time. I do not think that there is any workaround for transferring them to the HPC `scratch` as a desktop computer is simply not powerful enough to handle these amounts of data. For transferring data from A to B I use `rsync` which will keep the transferred file hidden until the transfer has been finished successfully. I prefer the following command:

rsync -axv --numeric-ids --progress -e "ssh -T -o Compression=no -x" *.bam user@path_to_hpc(...):/scratch/your_username/folder...

This will take time but `--progress` will give you a rough estimation for each file. Check if it might not be faster to directly download to the HPC. Talk to the people involved. Depending on how fast that drive is where the data are currently stored, a new download directly to `scratch` might save you quite some time.
Question: Hello everyone,

I'm working on improving the gene coverage for the chromosome of an E. coli strain, using existing WGS data. There are two WGS projects for this strain on NCBI (according to the assembly report, the assembly level for one of them is "scaffold", and for the other is "contig").

So far, I've obtained the contigs from both of these projects from NCBI (which I saved as two separate multifasta files), found a decent reference genome, and re-ordered the contig sets against it using Mauve. I then ran the re-ordered contigs through GLIMMER-3, which produced annotated Genbank format files for each WGS project contig set. Basically, I've been following the steps in this tutorial: http://www.microbialinformaticsj.com/content/3/1/2.

After plotting the two sets of annotated contigs against each other using DNA plotter, it looks like together they can give decent gene coverage. There are gaps in the CDS in each case which are complemented by the other.

I have been using contigs from both WGS projects so far, but I think I'm going to start using the scaffolds from the scaffold level assembly; as I've just found out that the scaffold level project includes two plasmids (the contig level WGS project is just a chromosome), and I know which scaffolds are the plasmids, so I can exclude them.

How can I go about assembling a single chromosome using data from both WGS projects? Is this possible? If not, which programmes can highlight the differences between the assembled chromosomes?

Should I be using the scaffolds instead of the contigs for the WGS project which is at a scaffold level?

Are there any programmes for closing gaps between/in scaffolds using additional contigs?

Thanks for your time,

Ronan

Answer: If the raw sequencing reads were also released, you could try GAM-NGS. Another option is to use MIRA to perform a genome-guided assembly, use the scaffolds as reference, and the two sets of contigs as two additional strains to be assembled. A third option is to use CAP3 to merge the assemblies, but I remember reading (from more than one source, though I do not recall a single one) this is not recommended, as it introduces mis-assemblies.

edit: I saw this post listed as similar to yours, seems like a good alternative to try.
Question: Dear Community,

scRNA-seq analysis has its own conventions. For example, CD56 is a canonical protein marker for the identification of Natural Killer (NK) cells. However, in scRNA-seq analysis, NCAM1 (gene for CD56) is not commonly used. Instead, GNLY and NKG7 are popular markers for identification of NK cell cluster in t-SNE or UMAP.

To address this issue, I have investigated several scRNA-seq datasets from different technologies such as 10X Chromium (v1, v2, v3) and Smart-seq2. It turns out that NCAM1 expression is very low in even FACS-gated (as CD56+) NK cells while GNLY and NKG7 are detected very well in this group.

My question is that why do some genes are more prone to dropout effect? What are the biological or technical explanations behind this issue?

Thank you for your contributions in advance.

Answer: The biggest factor is the amount of mature transcripts in a given cell. There may be some genes that are very quickly translated and/or encode very stable proteins, which would lead to relatively low numbers of free mature transcripts despite relatively high protein levels. Some genes are also transcribed in short bursts whereas others are more continuously transcribed, which would increase their chance of being present in the majority of cells at a given time point.

GC content may also play a role (very high or very low GC content will negatively impact PCR efficiency).

I can also imagine that some transcripts might be less amenable to poly-A-based capture, maybe due to secondary structures or additional interacting factors.

EDIT: This link provides a good run-down of the technical aspects limiting gene capture rates in general:

In short, the main reasons for drop-outs of individual genes in individual cells are due to (a) transcript abundance, (b) capture inefficiency, (c) amplification/sequencing bias.
Question: As we all know， bam file required by GATK need header in it. We can explicitly add it via the -R parameter.But，if you didn't use -R ,you also can use samtools -T genome.fa to add header. What is the difference between the two?

Answer: OK, now we can talk about this ;)

I would recommend you also have a look at the [sam specs][1] to understand what type of header informations (page 3ff.) exists in a bam/sam file.

The `-R` parameter during alignment adds a ReadGroup to the bam file. Doing this you can store informations like SampleName or used platform into the bam file. This is useful and often neccessary for downstream analyses like variant calling. It is best practice to always include a ReadGroup with at least the sample name.

What every aligner should automaticly store in the bam header is an informationen about the contigs (name and length) used during alignment. These information are introduced by `@SQ`. But this is not mandatory. But just like the ReadGroup there are tools that need this information. For example `samtools view`. If this information is missing you can use the `-T` parameter to provide the reference genome file used for alignment and `samtools` will extract the information it needs from there instead using the header.

fin swimmer

[1]: https://samtools.github.io/hts-specs/SAMv1.pdf
Question: Hi all,

I have a list of proteins that belong to an orthogroup, and id like to expand this gene family by searching for more members iteratively. I have read that jackhmmer in the HMMER package can do this however I am confused on the input. If i create an HMM profile and use that as input I get an error message:

> Error: Sequence file ortho1.hmm is empty or misformatted

So I am thinking of just giving it the protein sequences (not in a HMM profile) however I would need to ensure that it can make a profile with those sequences and use that to search iteratively, and output the final hmm model at the end. Anyone have experience with this and can give any advice?

Thank you!

Answer: It seems that you are feeding an HMM file to `jackhmmer`, which won't work. The error message tells you so. However, the alignment you used to build that HMM should be a valid input as long as it is in one of common formats (aligned FASTA, clustal, stockholm, etc).
Question: Hi guys,

I wonder if there're any up-to date publicly available genomes of bacteria, I could use as an input for Blast/Blat (fasta/2bit format) search. If so, would you send me links? I'm aware of database '16S ribosomal RNA sequences (Bacteria and Archaea)', however it's coverage is rather low. I appreciate your experiences.

EDIT: thank you guys, you've helped me a lot. It's a shame, that as an answer just one can be picked :)

Answer: I have a good blog post on this topic:

http://thegenomefactory.blogspot.com.au/2012/07/navigating-microbial-genomes-on-ncbi.html
Question: I have done genome assembly on a bacteria sample using SPAdes and Velvet. I have a contigs file from each of the assemblers. How can I compare these two files to see if their results are different or not?

Answer: I ended up using progressiveMauve, MUMmer (specifically MUMmerplots), QUAST, and clc_sequence_info (from CLC Genomics Workbench) altogether and I could come to some conclusions based on the results of all these methods.
Question: Im using two condition wild type and treated I get a scatter but i want to have two different colors Im not sure how to do it , I tried to group them by melting i got aesthetics error so i removed the gene column only i have two column but the data points are labelled as same

df <- HL_60
#gene = df[,1]
WT = df[,2:3]
ATRA = df[,4:5]
wt.mean = apply(WT,1,mean)
atra.mean = apply(ATRA,1,mean)
WT_ATRA = cbind.data.frame(wt.mean,atra.mean)
head(WT_ATRA)

library(ggplot2)
library(ggpmisc)
D1 <-WT_ATRA
head(D1)
my.formula <- y ~ x

p <- ggplot(data = WT_ATRA, aes(x= wt.mean,y = atra.mean) )+ geom_point(color = 'red',size = .9)+
geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
stat_poly_eq(formula = my.formula,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE)

The figure i get is this ![one][1]

Any suggestion or help how to put two colour for two different condition would be highly appreciated

[1]: https://i.imgur.com/AVdWuZK.jpg

Answer: The plot is in a 2D plane which means each point has two values associated with it (x,y) in this plane. in your case, each point is a gene and has (wt_mean,atra_mean) as coordinate. So, I don't get the point of coloring on the basis of wt_mean and atra_mean. What are you trying to infer?

The plot you mentioned in the comment using fold change as a factor of color code which can be achieved with the following code in `ggplot2`.

geom_point(aes(colour = cut(log2fc, c(-5, -1, 1, 5))),size = 2)+
scale_color_manual(name = "Fold Change",
values = c("(1, 5]" = "#FF6666",
"(-1,1]" = "#F5F5F5",
"(-5,-1]" = "#60B9FF"),
labels = c("<-2", "-2< & <2", ">2"))
`log2fc` is the log2 converted foldchnage information which is divied into three different bins (-5,-1), (-1,1),(1,5). Change them according to your rquirements.
Question: I am looking for pre-indexed BWA indices for human reference. Found a similar post with only a little help (some links are broken)

https://www.biostars.org/p/5474/

Answer: Illumina has a collection of "Ready-To-Use" reference sequences for commonly studied organisms which include the fasta reference sequence, Bowtie2 index, and BWA index files.

FIles: http://sapac.support.illumina.com/sequencing/sequencing_software/igenome.html?langsel=/in/

Documentation: ftp://ussd-ftp.illumina.com/README.txt
Question: Hi,

I have a vcf file with the following fields:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SampleX

The issue is that the ID column has '.' and no rsids. I have an annotation file in gtf format as follows:

1 AUGUSTUS CDS 106955274 106955412 0.89 + 0 transcript_id "g53770.t1"; gene_id "g53770";
1 AUGUSTUS intron 106955413 107007218 0.89 + . transcript_id "g53770.t1"; gene_id "g53770";

Does anyone know how I can use the annotation file to add rsids to my vcf file?

Thanks

Answer: you need a dbsnp-like vcf file, not a gtf

bcftools annotate -a dbsnp vcf.gz -c ID in.gz -O z -o output.vcf.gz
Question: Hi,

I may have a silly question. How should I get normalized counts from DESeq2? If I understand it right, rlog and vst are both normalized counts which one is better? or Should I just use "counts(dds,normalized=TRUE)" instead of rlog or vst. Thank you.

-X

Answer: For things like calculating fold-changes, one would normally use the output of count(dds, normalized=T). The other two you mentioned and so much count but rather transformed counts. These are useful for things like making a heatmap, or PCA, or anything that involves clustering/imaging.

BTW, none of these will give you the same fold-change that DESeq2 does, because you're not using a prior distribution.
Question: Hi All,

I have some PacBio data in bam format obtained from the Sequel platform.

I know that the PacBio sequencing is able to captures the methylation information, but I don't know how to "obtain" these information.

Is there a way to know which bases in the bam files are methylated using samtools?

Or is there a compleatly different way?

Thank You

Answer: I have complied some of the steps to do this in a previous post. I am not an expert on this by any means but most of PacBio's official documentation for this type of pipeline is outlined in the post linked below:

https://www.biostars.org/p/302290/#302647
Question: Hi, I am trying to train a predictor from some docking results. To do this I will need both true negative and true positive datasets. These datasets must have a PDB, small molecule and an indication of whether they interact.

I can probably make a true positive dataset by finding all PDBs which are bound to a small molecules but I'm not sure about the true negatives. I'm hoping that these datasets exist already from other publications.

Does anyone know where I might find these datasets? Thanks, Kevin

Answer: Most common method is to use DUD or newer version - DUDe.
The generated compounds - called decoys, active compounds can be found in ChEMBL.
Question: Hi all,

I'm trying to run HHblits but running into an error. I have several multiple sequence alignments that I have trimmed and concatenated and then converted the resulting .phy file to a .fasta (so HHblits can accept it).

Downloaded the database using:

wget http://www.gwdg.de/~compbiol/uniclust/2018_08/uniclust30_2018_08_hhsuite.tar.gz
tar xzvf uniclust30_2018_08_hhsuite.tar.gz

When I run HHblits though, I keep getting this error message:

$ hhblits -i sample.fasta -o hhblits_results.hhr -d ./databases/uniclust30_2018_08

- 13:51:18.331 ERROR: In /dev/shm/ebuser/avx2/HHsuite/3.0-beta.3/iccifort-2018.3/hh-suite-3.0-beta.3/src/hhdatabase.cpp:35: FFindexDatabase:
- 13:51:18.331 ERROR:	could not open file './databases/uniclust30_2018_08_cs219.ffdata'

my ./databases/uniclust30_2018_08 contains the following files:

uniclust30_2018_08_a3m_db	uniclust30_2018_08.cs219	uniclust30_2018_08_hhm_db	uniclust30_2018_08_md5sum
uniclust30_2018_08_a3m_db.index uniclust30_2018_08_cs219.ffdata uniclust30_2018_08_hhm_db.index
uniclust30_2018_08_a3m.ffdata	uniclust30_2018_08_cs219.ffindex uniclust30_2018_08_hhm.ffdata
uniclust30_2018_08_a3m.ffindex	uniclust30_2018_08.cs219.sizes uniclust30_2018_08_hhm.ffindex

Thank you

Answer: In the older versions at least, you had to provide the base name of the file itself, not just the directory.

With `hhsearch` its sufficient to use the `/path/to/databasename/databasename.ffdata`. With `hhblits` I believe it needs more than a single file, so you provide the basename of the files.

In your case it may look something like `./databases/uniclust30_2018_08/uniclust30_2018_08`.

This is also assuming you’re in the correct directory, since you’re using relative filepaths.
Question: Hello Biostars

I know that I can run blast with with "edible" output by using this command:

blastp -query -FASTA.faa - db MYDATABASE -evalue 1e-5 -outfmt "6 sseqid pident evalue bitscore staxids"

If I want to run the same command in parallel, i am using

ls *.fna| parallel -a - blastp -query {} -db MYDATABASE -evalue 1e-5 -outfmt 6 -out {.}.txt

My issue is when I am trying to combine the `outfmt` in parallel I got an error and I want to know if I can combine those two commands in one?

I tried:

ls *.fna| parallel -a - blastp -query {} -db MYDATABASE -evalue 1e-5 -outfmt "6 sseqid pident evalue bitscore staxids sscinames sskingdoms stitle" -out {.}.txt

This gives me:

Error: Too many positional arguments (1), the offending value: evalue
Error: (CArgException::eSynopsis) Too many positional arguments (1), the offending value: evalue

sorry for my terminology, but I'm a new kid on the block

Thank you in advance

Answer: I second the other suggestions of conflicts with how parallel and bash interpret quotations. I generally find it best to write a function for the command and then feed the function into parallel, avoiding exactly what you experience. I do not see a reason to use `-a` so leaving it out here (untested, based on your code, I do not know the blast syntax):

function blasty {

blastp \
-query "${1}" \
-out "${1%.fna}".txt \
-db MYDATABASE -evalue 1e-5 -outfmt "6 sseqid pident evalue bitscore staxids sscinames sskingdoms stitle"

}; export -f blasty

ls *.fna | parallel blasty {}

`MYDATABASE` is worth checking as well, is this a variable or in this case simply a spaceholder for the actual path to the file (or whatever it is, as said not a blast user myself)?
Question: Hello, everyone!

I have two types of files ".seq" and ".param". One has sequence info (obviously) and the other has parameters for PRIMER3.

I have about 100 files with those extensions.

For example,

```
File1.seq File1.param
File2.seq File2.param
File3.seq File3.param
File4.seq File4.param
File5.seq File5.param
```

I was trying to concatenate files with the same name using 'for loop' and 'cat', but it didn't seem to work. The code I was going to try

```
for I in ./*.seq ./*.param
do cat $i # I'm stuck here. I don't think $i can take both files.
```

Is there a way to load two variables in a single loop, or am I doing completely wrong?

Thank you!

Answer: One way is to get basename for one file and then append extension, e.g.

for F in *.seq, do N=$(basename "$F" .seq); cat $F "$N".param > "$N".out; done