cannin/biostars_qa huggingface dataset:
Question: Hi,

I was wondering whether it's possible to have a heatmap in R with two different colour scales.

I have two experiments with two different results of a different scale. The ID list list is the same. I would like to plot them both on the same image if possible. But as they are of different values, I can't plot them on the same scale.

I would like the heatmap to have a different scale and a different colour-space. The fist one is from -10 to 10 and should have e.g. the colour scale green-red. the second is from -3 to 3 and should be plotted in yellow-blue colour space.

Is something like that possible?

Thanks

Assa

Answer: Yes it is possible to plot using heatmap or heatmap.2 with a different color scheme for each column, but not without some trickery.

First, you have to scale your columns to unit variance, because it would be a mistake not scale data used for clustering if we already know they are on a different scale. Most distance measures are sensitive to differences in variance between the variables and your first experiment has a much larger range than the second one. If both distributions are approximately normal or at least central, it is reasonable to assume that the first variable has also a larger variance and will therefore dominate the cluster analysis.

library(gplots)
### make some example data
x = matrix(c(runif(10, -10, 10), runif(10,-3,3)), ncol=2)
### scale the matrix
xs = scale(x)
### save the dendrogram:
row.dend = as.dendrogram(hclust(dist(xs)))
### make a single color scale red-black-green, yellow-white-blue
mycols = c(colorRampPalette(c("red","black","green"))(10), colorRampPalette(c("yellow","white","blue"))(10))
# did I say that this was ugly?
# adjust columns of xs to plot:
xs[,1] = xs[,1] - max(xs[,1]) - 0.1 # make sure they are far away enough from 0!
xs[,2] = xs[,2] - min(xs[,2]) + 0.1
## control color breaks are centered at 0 and override automatic color breaks,
## which can result in use of wrong color on either side depending on distribution
breaks = (seq (-1,1, length=21))
## plot:
heatmap.2(xs, col = mycols, Rowv = row.dend, dendrogram="row", scale="none", breaks=breaks)
## compare this to the outcome using a single color scheme which is much more intuitive
heatmap.2(x, dendrogram="row", scale="column", col=colorRampPalette(c("yellow","white","blue"))(20))

Oh, and did I say that this solution stinks (you have to adjust color generation, or generate color breaks manually), even more so if you try to extend it to more than two colors. You get all sorts of problems with the values close to 0 now, also the cluster result doesn't look very convincing now, but in the end, this plot shows, why you shouldn't use different color-schemes in a single heatmap. (Also there is a difference in the dendrogram between both approaches, need some time to find out why).

The best way would be to make a derived version of `heatmap.2` or `image`, that supports a list for its `col` argument, allowing to pass one color vector for each column, please go ahead ;P

**tl;dr: Yes it is possible, but you don't want to do it.**
Question: If you look at the materials and methods of ([https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5449402/][1]) this paper, they have suggested using any three tools from each method (sequence homology-based methods, supervised learning methods, protein sequence and structure-based methods and consensus-based methods) to find the deleterious SNP. How do I choose the three best tools for each method?

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5449402/

Answer: Take your pick...

#Sequence and evolutionary- / homology- / conservation-based methods (for missense predictions)

- SIFT
- Align-GVGD

- Mutation Assessor
- PANTHER
- MAPP
- MUpro
- MutationTasster
- I-Mutant
- KD4v
- FATHMM
- dbNSFP

#Protein sequence and structure-based methods (for missense predictions)

- PolyPhen-2
- LS-SNP/PDB
- SNPeffect 4.0

#Supervised-learning methods (for missense predictions)

- PMut
- SNAP
- PhD-SNP
- SNPs&GO
- Parepro
- CanPredict
- nsSNPAnalyzer
- MutPred2
- HANSA
- MutationTaster

#Splice predictions

- NNSplice
- GeneSplicer
- MaxEntScan
- HSF
- NetGene2

#Protein modelling (from amino acid sequence)

- Protein
Model Portal

*[uses various modelling algorithms and produces PDB files, which can be loaded into protein viewers like Jmol]*

#Non-coding (i.e. regulatory)

- CADD (germline
variants)
- DANN
(germline variants)
- FATHMM-MKL
(germline variants)
- GWAVA
(germline variants | somatic mutations)
- Funseq2 (somatic
mutations)
- SurfR
(rare variants | complex disease variants | all other variants)

# Other

- GWAS3D / GWAS4D

-------------------------

For further reading:

- UK perspective, refer to the guidelines by the Association for
Clinical Genomic Science:
http://www.acgs.uk.com/quality/best-practice-guidelines/
- US perspective, refer to American College of Medical Genetics:
https://www.acmg.net/docs/Standards_Guidelines_for_the_Interpretation_of_Sequence_Variants.pdf

------------------------------------------------------------
------------------------------------------

I assure you that there are even more tools than these. These are categorised to the best of my and my colleague's knowledge.

Previous posting: https://www.biostars.org/p/286364/#286483

Kevin
Question: $ awk -vPRIMER="ACTG" '{if ($0~/^>/) { print $0; } else { print PRIMER$0; }}' in.fa > out.fa

(by Alex Reynolds)

Can this code be modified to add primer sequence in the end of the complete sequence, what we have generated in the previous code? And I assume we will already use the reverse complement sequence so that we dont have to include that thing in the code.

Answer: Reposition `PRIMER$0` to `$0PRIMER` to move the string to the end.
Question: I have downloaded a histone mark data from roadmap epigenomics:
CD4+_CD25-_CD45RO+_memory_primary_cells/H3K4me3/GSM772862_BI.CD4+_CD25-_CD45RO+_Memory_Primary_Cells.H3K4me3.Donor_62.bed

but I am not quite sure I understand the content. It is chip-seq analysis, with peak calling, so I expected to see coordinates of the peaks and their height.

However I see such information in a file: chromosome, start, end, id, something, strand

```
chr1 9984 10183 B09JPABXX110526:5:1101:19182:47634 0 -
chr1 9990 10189 B09JPABXX110526:5:2207:9781:41112 0 -
```

I was intended to use this data to generate plots with GVIZ to get such plots but I am not sure if start and end positions will be enough somehow.

Answer: Ah, Roadmap data. It's a bit of a head-ache to understand it sometimes, but bear with me:

What you are looking at is raw (unconsolidated) input that was output by Pash mapper. It is a misformatted [bed file][1]. All values in start column should have one subtracted from their coordinate (not respecting the strand), i.e. it should be:

```
chr1 9983 10183 B09JPABXX110526:5:1101:19182:47634 0 -
chr1 9989 10189 B09JPABXX110526:5:2207:9781:41112 0 -
```

I am not really sure about what name parameter encodes or why the score is zero...

Once you fix this (using [bedtools slop][2], `bedtools slop -l 1 -r 0 -g hg19.genome`), just pass the fixed bed-file to some sort of pileup tool (for instance, MACS).

The way ROADMAP does it:

1. They shorten all the reads back to 36 (so its consistent across all experiments)
2. They filter duplicate reads, and [reads that could not be mapped][3] to the genome at all had they been 36bp long
3. They estimate fragment length using [SPP][4] and run [macs2][5] with appropriate fragment lengths (`--nomodel --ext-size=fragment_length`) to generate both peak lists and the signal/foldchange track.

I really suggest you skip doing the preprocessing yourself and just get the data from their [download page][6]

Namely you want to look at c (peak calling) or section d (signal tracks). Either go with consolidated reads for your cell line (consolidated = all technical/biological replicates lumped together, recommended), or with unconsolidated (which is what you are looking at at the moment).

[1]: https://genome.ucsc.edu/FAQ/FAQformat.html#format1
[2]: https://bedtools.readthedocs.org/en/latest/content/tools/slop.html
[3]: https://sites.google.com/site/anshulkundaje/projects/mappability
[4]: https://code.google.com/p/phantompeakqualtools/
[5]: https://github.com/taoliu/MACS/
[6]: http://egg2.wustl.edu/roadmap/web_portal/processed_data.html
Question: Hello,

I understand that an important aspect of eQTL analysis is accounting for confounding variation in the expression data, as many factors can affect gene expression, including sex, age, lifestyle factors etc. Including these factors as covariates in the model therefore increases the chances of identifying genuine eQTL effects.

I have RNASeq data for approximately 100 individuals, along with corresponding genotype data for approx 600,000 SNPs. I'd like to use MatrixeQTL to do an eQTL analysis using this data. I'm also in the fortunate position of having a lot of clinical data for these subjects, describing a variety of things, including the usual suspects: age, gender, ethnicity etc. All in all the clinical data comprises 200 fields. However, as is usually the case with clinical data, much of it is incomplete (i.e. contains NAs).

Any number of these fields could be influencing gene expression, but i'm unsure how best to account for them in the analysis. Do I include all them as covariates? Only the "obvious" ones (age, sex etc)? Or is a better approach to use principle components of the expression data as covariates?

Any guidance would be appreciated.

Answer: Dont think there is a golden standard approach here, much will depend on what you research question is?
Couple things you can do:

1) impute your clinical data to fill in the NA

2) make PC's from your RNAseq data and check for correlation between these PC's and any of your clinical covariates

3) some people prefer to unbiasly normalize their RNAseq data using PEER (or similar) software to adjust for hidden covariates as well.
Question: Hello BioStar,

After some time working with Illumina and pipelines, I've identified a bottleneck when getting early "draft" results from a run. Figuring out how the sequencing data looks like in real time, as opposed to wait for the run to finish after ~11 days.

The goal would be to get an estimate of how many reads one could expect to get for each sample, thereby guiding setup for a subsequent run for topping up the data for those samples that do not reach the required amounts. It would be help a lot if we could reduce the the wall-clock time for reaching a decision on which samples need to be re-run.

When it comes to implementation, I've been thinking on a file status daemon such as Guard[1], coupled with the CASAVA/OLB tools from illumina, performing basecalling and demultiplexing as soon as the files get written to disk, without having to wait for the whole run to finish.

Other more high-tech solutions would involve Hadoop Seal[2] and Flume[3] to do the streaming part into Hadoop... but before digging more into this issue I wondered what are you guys doing to get a draft view on a running sequencing run.

Cheers & happy new year Bio* !

[1] https://github.com/guard/guard

[2] http://biodoop-seal.sourceforge.net/installation_generic.html

[3] http://www.slideshare.net/cloudera/inside-flume

Answer: We're waiting for the whole run to finish before getting any QC information. Optimizing that step like you say would surely be beneficial, and if you can get the base calling on the index read to run as soon as it's been sequenced, then I think the rest would be simple.

Since you only want to get an estimate of the read counts for each of the multiplexed sample, I would try simply using a script to calculate the index read frequencies. You could even exclude reads that have quality problems from your counts. I don't think performance would be an issue since you'd only be dealing with the index reads for one run but, if necessary, you could quickly copy the index reads to a Hadoop cluster and run a custom Pydoop script (or even recycle the Hadoop workcount example) to generate a table (index read, count), which you could join to your sample sheet to get your sample read counts.

HTH.
Question: I am using the PRSice tool to calculate PRS scores for a dataset (d1) using a set of weights from the PGS catalog. Although the weights were generated from a different cohort (d2), it appears that the original GWAS was based on the same dataset (d1) that I'm using. Does this method inflate the predictive qualities of the PRS for the cohort in d1? Or should we even use it in that case? (Still very new to the concept of PRS scores)

Answer: Unfortunately, if the same data set is used for the GWAS estimate, you can't include that in the PRS analysis, or you will get a highly inflated R2 and significance that cannot be replicate in other data set.

Do you mean that d1 is part of d2? If d1 is a subset of d2, one method is to use the inverse meta-analysis approach. you should be able to find the equation here: https://www.biorxiv.org/content/10.1101/203257v1.full.pdf
Question: The Bam format was developed at a time when I/O (disk / network) was not a limiting factor for genomics data analysis. Currently it is. Especially when you have pipelines like GATK that process BAM files in multiple steps and read and write bam files for every step (instead of in memory streaming).

I know raw sequencing data (SOLiD & PacBio) and other large (scientific, financial etc.) datasets are stored in [HDF5](http://www.hdfgroup.org/HDF5/). And there are initiatives like [CRAM](http://www.ebi.ac.uk/ena/about/cram_toolkit) and [reduced bam](https://www.broadinstitute.org/gatk//events/2038/GATKwh0-BP-4-Compression.pdf).

Still I believe most organizations are using regular bam files. Why is this? Is there a roadmap towards development (and adoption!) of improved formats for storing DNA alignment data?

Answer: > Why is this?

because, most tools (samtool 0.1.19), workflows... still use BAM.

you never change a winning team.

HDF5: I don't think it's the best tool to store data (you want Streaming , you want to avoid padding of data if the record is tool small, you want bgzip compression,... ). On a personnal side, I think HDF5 is a nightmare to program in C. https://www.biostars.org/p/69/

oh and why don't people use XML or JSON instead of VCF? why don't people all use the same best programming language (which is java *) :-)

(*) troll.
Question: What's ab initio method in bioinformatics? Can you give me a general overview of this? You can indicate me a source where can find this informations.

Answer: As https://www.biostars.org/u/24526/ said, there are several different "*ab initio*" in bioinformatics, e.g.:

- *ab initio* protein folding
- *ab initio* gene prediction
- *ab initio* genome assembly

And so on. In general, *ab initio* in bioinformatics means "without direct evidence". For example, *ab initio* genome assembly is assembling a genome from sequencing reads alone, without any other information like a closely related genome or a genetic map.
Question: I have about 30 GB (gzipped fastq) of small RNA (miRNA, etc.) sequencing data that I've been asked to align to the human genome (hg38), reporting for each read *every* position in the genome where it matches exactly. Bowtie2 has been my go-to mapper for short reads in the past, so I put together a pipeline something like this:

zcat Sample_1.fastq.gz | \
cutadapt -a TGGAATTCTCGGGTGCCAAGG -f fastq --overlap 15 \
--trimmed-only --minimum-length 10 - | \
bowtie2 -x ~/references/hg38/hg38-bt2/hg38 -U - -q --phred33 \
--end-to-end --very-sensitive -a --score-min C,0,-1 -L 10 \
picard-tools SortSam INPUT=/dev/stdin OUTPUT=Sample_1.bam;
picard-tools BuildBamIndex INPUT=Sample_1.bam OUTPUT=Sample_1.bam.bai

However, the Bowtie2 manual carries [this warning][1] about the `-a` option for reporting all alignments for each read:

> Some tools are designed with this reporting mode in mind. Bowtie 2 is not! For very large genomes, this mode is very slow.

Nevertheless, because I already had it set up, I decided to give it a try. Apparently, the statement is quite accurate, because one sample (about 300 MB before trimming) has already been running for 168 CPU-hours without finishing, which is at least an order of magnitude or two longer than I'd typically expect a 300 MB fastq file to require when looking for just one or a few alignments for each read. So, is there either another set of options for bowtie2 or, more likely, another aligner better suited for the task of quickly finding all exact matches for each read?

[1]: http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#a-mode-search-for-and-report-all-alignments

Answer: Just to add, if all you care about is all perfect alignments, I recommend these BBMap flags:

**ambig=all vslow perfectmode maxsites=1000**

It should be very fast in that mode (despite the vslow flag). Vslow mainly removes masking of low-complexity repetitive kmers, which is not usually a problem but can be with extremely short sequences like microRNAs.
Question: Hi,
I have a list of partial seq IDs, and a fasta file, and I want to get a fasta file that contains all the sequences that match the IDs in the list.

Example:

list_ID.txt

> AF306937.1
> S56354.1

Requested output:
proteins_output.fasta

> lcl|AF306937.1_prot_AAG42454.1_1
> MGMDPIALQAGFDLLGDGRPETLWLGIGTLLMIIGTFYFIAQGWGVTDKEAREYYAITILVPGIASAAYLAMFFGIGVTEVELASGAVLDIYYARYADWLFTTPLLLLDLALLAKVDRVSIGTLIGVDALMIVTGLIGALSKTPLARYTWWLFSTIAFLFVLYYLLTSLRSAAAQRSEEVQSTFNTLTALVAVLWTAYPILWIVGTEGAGVVGLGVETLAFMVLDVTAKVGFGFALLRSRAILGETEAPEPSAGT
>
> lcl|S56354.1_prot_AAB19870.2_1
> MDPIALQAGFDLLNDGRPETLWLGIGTLLMLIGTFYFIARGWGVTDKEAREYYAITILVPGIASAAYLAMFFGIGVTEVELASGTVLDIYYARYADWLFTTPLLLLDLALLAKVDRVTIGTLIGVDALMIVTGLIGALSKTPLARYTWWLFSTIAFLFVLYYLLTSLRSAAAKRSEEVRSTFNTLTALVAVLWTAYPILWIVGTEGAGVVGLGIETLAFMVLDVTAKVGFGFVLLRSRAILGETEAPEPSAGADASAAD

I wrote this loop for bioawk, but it isn't quite working as I want.

When I run this loop, I get an empty "proteins_outout.fasta" file:

while read p; do
bioawk -c fastx '$name ~ /$p/ { print ">"$name"\n"$seq"\n"; }' proteins_db.fasta > proteins_outout.fasta;
done "$name"\n"$seq"\n"; }'
proteins_db.fasta

>lcl|AF306937.1_prot_AAG42454.1_1 MGMDPIALQAGFDLLGDGRPETLWLGIGTLLMIIGTFYFIAQGWGVTDKEAREYYAITILVPGIASAAYLAMFFGIGVTEVELASGAVLDIYYARYADWLFTTPLLLLDLALLAKVDRVSIGTLIGVDALMIVTGLIGALSKTPLARYTWWLFSTIAFLFVLYYLLTSLRSAAAQRSEEVQSTFNTLTALVAVLWTAYPILWIVGTEGAGVVGLGVETLAFMVLDVTAKVGFGFALLRSRAILGETEAPEPSAGT

Why is my loop not working as I want?

Would appreciate any help!

Thanks

Answer: I think a) you have to pass the shell variable `$p` to `bioawk` properly using `-v` and b) you are currently overwriting your output file in every iteration as you use `>` rather than `>>`:

while read p
do
bioawk -v var="${p}" -c fastx '$name ~ var {print ">"$name"\n"$seq}' proteins_db.fasta >> proteins_outout.fasta
done < list_ID.txt
Question: I have sequencing data of five bacteria, which were generated using Illumina MiSeq. Four of them were sequenced using a paired-end 2x300 protocol and one was sequenced using the nextera mate-pair protocol.

My question is: What are the softwares that you recommend me to assemble these genomes (the largest has almost 8 Mbp)?

I have access to a CLC Workbench. It seems quite ease to use, but I dont know if it is the best one. Most of papers that I found that evaluate the performance of assemblers from two-three years ago.

I also have to mention that I have two i7 with 8 threads PCs available for this objective (one with 32 and another with 8 Gb RAM).

Thanks in advance.

Answer: For bacterial genome, [GAGE-B paper][1] (2013) compare 8 genome assemblers:

- ABySS v1.3.4
- CABOG v7.0
- MIRA v3.4.0
- MSRCA v1.8.3
- SGA v0.9.34
- SOAPdenovo2 v2.04 + GapCloser v1.12
- SPAdes v2.3.0
- Velvet v1.2.08

All GAGE-B data and assembly recipe are available [here][2].

For more recent comparison of genome assemblers, have a look [here][3].

As each bacteria genome size and GC% is different, you need to check [these reproducible Benchmarks][4].

[1]: http://bioinformatics.oxfordjournals.org/content/29/14/1718
[2]: http://ccb.jhu.edu/gage_b/index.html
[3]: http://nucleotid.es/
[4]: http://nucleotid.es/benchmarks/
Question: Dear Community,

I am trying to use bowtie2 or bwa multi-core/threads option to align my reads to genome (-p 4). I am using a MacBook pro m2 chip and it seems cannot fully parallelly use all the cores. using 'top", it shows 4/1 for threads (#TH) and only one process for bowtie2.
Could someone clarify this? and let me know how to enable parallel multi-threads for using these aligners on Mac M2 chips?
Thank you!

-Xiaokuan

Answer: Programs do not use all the threads/cores all the time. They will do so during appropriate part of the analysis. As long as you have enabled use of multiple threads on your command line there is nothing else you need to do.
Question: Dear Biostars,

I have a NGS data for thousands of people merged into a multisample vcf which has been split into genomic chunks (roughly 1000 chunks). These have been annotated with VEP. One of those annotations is the LOFTEE "high confidence" marker.

I want to count per sample the number of LOFTEE "HC" hits there are per sample across all chunks so the final output would simple be the sample ID and a number representing the total number of high confidence markers that person has.

I have tried:

bcftools view -S $sample_list | \
bcftools +split-vep \
-d \
-f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Feature\t%Consequence\t%Existing_variation\t%LoF\n' \
-i "LOF='HC'" \
-o output_${i}.txt

where $sample file is an array over all the chunks but I am not getting it per sample ID.

Any help would be greatly appreciated.

Answer: Correct me if I'm wrong, but the LOFTEE annotation depends on the variant and not on the genotype, therefore the `LOF='HC'` annotation would be inside the INFO column. If that'd the case, I'd suggest a search+split+count strategy rather than a split+search+count one.

Firstly extract `LOF='HC'` variants to a new multisample file

bcftools view -i ''LOF='HC'" -Oz -o multisample.LOF_HC.vcf.gz multisample.vcf.gz

and secondly split that new file by sample and count private variants in it

for sample in `bcftools query -l multisample.LOF_HC.vcf.gz`; do
echo -n "$sample "
bcftools view -c1 -s $sample multisample.LOF_HC.vcf.gz | grep -cv ^#
done > samples.txt

Note that this code is not tested, and that the split loop can take quite a while if the number of samples is large. If that'd be the case, then I'd go for a programmed solution by parsing the file directly:

zgrep -v ^## multisample.vcf.gz \
| perl -lane 'if (/^#/) {
foreach (@F) { push @h, $_; $s{$_} = 0 }
} elsif ($F[7] =~ /LOF='HC'/) {
for ($i = 9; $i <= $#F; $i++) {
($gt = $F[$i]) =~ s/:.+//;
if ($gt =~ /[1-9]/) { $s{$h[$i]}++ }
}
} END {
foreach (@h[9..$#h]) { print "$_\t$s{$_}" }
}' > samples.txt

(This code has been successfully tested on 1000g multisample files)
Question: I have a list of gene name in a file

```
CHRNB2
EGR2
GCK
KRT14
LMNA
FGF3
TK2
ABCC8
```

How can U map them to uniprot ID?

**P.S** I tried Uniprot "ID mapping" (from-"GENEID" to-"UNIPROTKB AC").But it couldn't map.

Please suggest me what to do.Thnx

Answer: These look to be gene symbols not [Entrez Gene][1] identifiers, so identifier mapping for Entrez Gene would not work. However UniProt include most gene symbols and their various synonyms in their data, so a query should work and find the set of UniProtKB entries which match. For example:

1. Go to the [UniProt.org website][2]
2. Select "Protein Knowledgebase (UniProtKB)" for "Search in"
3. For the "Query" enter 'gene:' followed by the gene symbol (e.g. "gene:CHRNB2")
4. Click the "Search" button
5. The results contain a list of UniProtKB entries matching this gene symbol (e.g. [for CHRNB2][3])

I suspect you may have a specific species in mind so you might want to use additional terms to limit the results further.

Once you have worked out the form of the required query you could use the [UniProt.org REST API][4] to script the required queries and return only selected data.

[1]: http://www.ncbi.nlm.nih.gov/gene/
[2]: http://www.uniprot.org/
[3]: http://www.uniprot.org/uniprot/?query=gene%3ACHRNB2
[4]: http://www.uniprot.org/faq/28