cannin/biostars_qa huggingface dataset:
Question: This is a tiny script that is enough to get gcc to fail when trying to build:

#include
#include
#include
#include

int main(int argc, char *argv[]){

samFile *fp_in = hts_open(argv[1],"r"); //open bam file
bam_hdr_t *bamHdr = sam_hdr_read(fp_in); //read header
bam1_t *aln = bam_init1(); //initialize an alignment
}

This is my minimal command for trying to build:

gcc test.c /mnt/work/me/software/anaconda/pkgs/htslib-1.9-hc238db4_4/lib/libhts.a -I /mnt/work/me/software/anaconda/pkgs/htslib-1.9-hc238db4_4/include/

What happens is:

/mnt/work/me/software/anaconda/pkgs/htslib-1.9-hc238db4_4/lib/libhts.a(hts.o): In function `decompress_peek':
hts.c:(.text+0x14a5): undefined reference to `inflateInit2_'
hts.c:(.text+0x14c6): undefined reference to `inflate'
hts.c:(.text+0x14f9): undefined reference to `inflateEnd'
/mnt/work/me/software/anaconda/pkgs/htslib-1.9-hc238db4_4/lib/libhts.a(cram_io.o): In function `itf8_decode_crc':
cram_io.c:(.text+0x18e0): undefined reference to `libdeflate_crc32'
cram_io.c:(.text+0x1938): undefined reference to `libdeflate_crc32'
cram_io.c:(.text+0x19b6): undefined reference to `libdeflate_crc32'
cram_io.c:(.text+0x1a5a): undefined reference to `libdeflate_crc32'
cram_io.c:(.text+0x1b2d): undefined reference to `libdeflate_crc32'

1) What am I missing to get this script to compile?

2) And how do I find out where htslib is installed on the end users system? I can am using setup.py to install the script so I am free to run arbitrary code to find it.

Answer: You are moving from building existing packaged tools, which is (admittedly less than ideally) something that many bioinformaticians need to do when they'd really like to be just using the tools, to writing the build infrastructure from scratch for your own tools, which is a software development task and in a different category.

So the HTSlib etc maintainers go to great lengths to make building HTSlib, SAMtools, etc Just Work, but there is less pressure to provide full instructions and Just-Workingness for people writing their own tools that use HTSlib, as there is a greater expectation that such people know what they're doing.

To get this **program** to compile, you are missing `-lz` from your link command. And something similar for the other errors. Unfortunately one of the skills you are going to have to learn in order to write a tool like this is the skill of debugging build issues like this one. There aren't many shortcuts other than experience, and help and advice from a more experienced programmer or sysadmin at your institute who can look over your shoulder at your screen.

The way to gain that experience starts with googling those error messages to figure out what you need to add to get this to compile and link. And learning about the different stages (compiling, linking, installing) involved in building a C program.
Question: Hi friends,

I am trying to download ht-seq file form GDC portal but it gives me one file for each patients. Can you please let me know how to download one file including all patients together for all 60000 genes?

Is there any R code for this?

Thank you

Answer: Whatever method you tried, internally it stills need to download all files, before the RNA API is in production.

There are some R tools doing that, or you can use "GDC RNA tool" listed here https://gdc.cancer.gov/access-data/gdc-community-tools for aggregation
Question: Hi all,

I would like to know if my approach makes sense or not.

I have a count of histone mark peaks (Encode BigWig files) in an 80Kb window from hg19.

I have sampled 1000 positions (80Kb in length) of hg19 and generated a null empirical distribution based on the counts of peaks from each of the 1000 positions. All looks good and we have a normal looking distribution.

Now to check if my original 80Kb window is significantly enriched I can simply check to see if my observation counts falls inside or outside of the 95% of the null distribution count data by looking at the mean and standard deviation of the null distribution (95%-99% rule). I would also like to get an approximate P-Value so this is where I am slightly unsure. I think I can do two of the following depending on a one tail or two tail test:

***One Tail Test***: Get the total number of observations from the null distribution that are greater than or equal to my observation and divide by 1000.

***Two Tail Test***: Subtract my observation count from the mean of null distribution count and then take the total number of observations greater than or equal to this absolute value difference and divide this by 1000.

Conformation or Suggestions will be greatly appreciated!

Answer: Yes, that's how it should be done. In both cases you're testing for the fraction more extreme than what you observed. In one case, you only care about "more extreme in one direction" (a 1-tailed test) while in the other you care about "more extreme in either direction" (a 2-tailed test).
Question: Hi,

I'm trying to export the log normalized count data from Seurat in to a .csv file and get the following error:

library (Seurat)
fib.data <- subset(onlyharmonfib, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mito <5)
fib.data_norm<- NormalizeData(fib.data)
fib_count <- GetAssayData(object = fib.data_norm, assay = "RNA", slot = "data")
write.csv(as.matrix(fib_count), file = "file_path", quote = FALSE)

Error in file(file, ifelse(append, "a", "w")) : cannot open the connection

Is there an easier method to get the data in a .csv file?

Thank you for your help!

Answer: I would use the `data.table` library to write the file, otherwise it could take a while to write to disk.

library("data.table")

fwrite(as.data.table(fib_count, keep.rownames="feature"), "counts.csv")
Question: Do you know any public scientific SQL server ?

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29

ENSEMBL http://uswest.ensembl.org/info/data/mysql.html

GO http://www.geneontology.org/GO.database.shtml#mirrors

(I'll give a +1 to each correct answer)

Answer: Ensembl and EnsemblGenomes both have MySQL access to the main dbs and mart dbs:

mysql -hensembldb.ensembl.org -P5306 -uanonymous
mysql -hmysql.ebi.ac.uk -P4157 -uanonymous

MySQL DB (BioMart)

mysql -hensembldb.ensembl.org -P5316 -uanonymous

EnsemblGenomes

mysql -hmysql.ebi.ac.uk -P4157 -uanonymous
Question: Dear all,

I was recently asked in a lab meeting whether there may be a gene length bias in the results produced by WGCNA, if the input data consists of normalized counts or variance-stabilized transformed data from DESeq2 (which does not correct for gene length). I am not sure how to answer this, as [there are several papers that used DESeq2 normalized counts as input for WGCNA][1], and I thought that this was actually [recommended by the authors of WGCNA][2].

In other words, could my detected expression modules be biased for gene length, with some modules being particularly driven by the length of the genes in that module as opposed to its actual expression?

[1]: https://scholar.google.com/scholar?hl=en&as_sdt=0%2C9&q=deseq2%20wgcna&btnG=
[2]: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html

Answer: Like DEG analysis, WGCNA is looking at each gene independently of one another. All that matters is how the expression of the same gene compares across samples. That comparison will be proportionally the same whether you normalize by gene length or not, because all the values within each comparison are from the same gene, so they have the same length. The modules are formed by finding genes that have the same patterns of expression across samples, not by comparing actual gene expression values between different genes.
Question: Hi

When using following code in R, I am getting error

lolliplot(mutation_data, fillCol="mutation_type", labelCol="amino_acid_change", proteinColour = "#ffffff", paletteB ="green", txtSize=3, txtAngle=20, host="dec2013.archive.ensembl.org")

Using the following host: dec2013.archive.ensembl.org for biomaRt queries to change the ensembl annotation version alter this parameter!
Querying biomaRt for transcript sequence
Request to BioMart web service failed.
The BioMart web service you're accessing may be down.
Check the following URL and see if this website is available:
http://dec2013.archive.ensembl.org:80/biomart/martservice?type=registry&requestid=biomaRt
Error in if (!grepl(x = registry, pattern = "^\n*")) { :
argument is of length zero

How to solve this?

Answer: That particular Ensembl archive is no longer available. You can view the list of available archives by following the link in the error message (you'll be redirected to an appropriate page) or by using the `listEnsemblArchives()` function in **biomaRt** e.g.

> biomaRt::listEnsemblArchives()
name date url version current_release
1 Ensembl GRCh37 Feb 2014 http://grch37.ensembl.org GRCh37
2 Ensembl 95 Jan 2019 http://jan2019.archive.ensembl.org 95 *
3 Ensembl 94 Oct 2018 http://oct2018.archive.ensembl.org 94
4 Ensembl 93 Jul 2018 http://jul2018.archive.ensembl.org 93
5 Ensembl 92 Apr 2018 http://apr2018.archive.ensembl.org 92
6 Ensembl 91 Dec 2017 http://dec2017.archive.ensembl.org 91
7 Ensembl 90 Aug 2017 http://aug2017.archive.ensembl.org 90
8 Ensembl 89 May 2017 http://may2017.archive.ensembl.org 89
9 Ensembl 88 Mar 2017 http://mar2017.archive.ensembl.org 88
10 Ensembl 87 Dec 2016 http://dec2016.archive.ensembl.org 87
11 Ensembl 86 Oct 2016 http://oct2016.archive.ensembl.org 86
12 Ensembl 85 Jul 2016 http://jul2016.archive.ensembl.org 85
13 Ensembl 84 Mar 2016 http://mar2016.archive.ensembl.org 84
14 Ensembl 83 Dec 2015 http://dec2015.archive.ensembl.org 83
15 Ensembl 82 Sep 2015 http://sep2015.archive.ensembl.org 82
16 Ensembl 81 Jul 2015 http://jul2015.archive.ensembl.org 81
17 Ensembl 80 May 2015 http://may2015.archive.ensembl.org 80
18 Ensembl 79 Mar 2015 http://mar2015.archive.ensembl.org 79
19 Ensembl 78 Dec 2014 http://dec2014.archive.ensembl.org 78
20 Ensembl 77 Oct 2014 http://oct2014.archive.ensembl.org 77
21 Ensembl 76 Aug 2014 http://aug2014.archive.ensembl.org 76
22 Ensembl 75 Feb 2014 http://feb2014.archive.ensembl.org 75
23 Ensembl 67 May 2012 http://may2012.archive.ensembl.org 67
24 Ensembl 54 May 2009 http://may2009.archive.ensembl.org 54
Question: I have a multigenome vcf file. Suppose the file has samples A to Z, but I want to extract the subset of samples B to G and extract a small vcf file. How can I make such subset vcf file?

Answer: from https://samtools.github.io/bcftools/bcftools.html#view

bcftools view -s samplelist

or

bcftools view -S samplefile

would do the job. docs are your friends ;)
Question: I am developing machine learning models to classify disease/non-disease patients using gene expression data. I have applied LASSO to select features and built classifiers using some of the top features after feature selection. Now I have to do external validation on an independent test set to judge my model's generalizability. The problem I am facing is while doing this part.
The training set is built using a GEO dataset that is built on [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform. However, the test set of my interest that has the same disease/non-disease labels as my training is an Array Express dataset built on A-MEXP-2210 - Illumina HumanHT-12_V4_0_R1_15002873_B platform. So, some of the top features that I would have selected from my test set to validate my model, is missing altogether. What should be the ideal way to validate here?

1. Use only those genes that are there in test set, select those genes from training set and build a model?
2. KNN-impute those genes in the test set and do the analysis?
3. Assign expression value zero for those genes in the test set and do the analysis?

Answer: I don't recommend imputing (zero, knn, or otherwise). If those features are missing altogether in the test set, you can't really infer what values those features (if they were to exist) would have unless you have some prior knowledge (which should not be taken from your training set) about what values those features should have.

The best way is to intersect the features from the training set and testing set, and work off of that for building your model.

If you're afraid of your model losing important features as a result, but still want to validate the model on an independent dataset, then I would go with imputing. Since entire features are missing, you wouldn't be able to use KNN (and you certainly shouldn't use any information at all from the training set for imputation purposes) -- I'd probably just assign the missing features to be the average expression of all genes for that sample. Just don't expect your model to perform as good since you're giving some semi-arbitrary values to some features.

Another way to do validation (if you're only interested in feature selection, i.e. if you just want to see whether you've selected "good" features that can predict disease vs. non-disease) is to take the top features you found, use only the ones that also exist in the second dataset, and then train + cross-validate your model using just the second dataset. This will not tell you whether the model you trained using the first dataset is generalizable to other datasets, but it will tell you whether your features (which you selected using machine learning on your first dataset) can be predictive of disease/non-disease in another dataset.
Question: Hey,
I have the following problem: using the `karyoploteR` package in R, i want to plot chromatin configuration along the genome as a colored bar. The GRanges Object with the chromatin regions has a column for that ("V9") with rgb values associated with each region. However, the `kpPlotRegions` function requires the colors in hex-format. I tried to export the Object, change rgb to hex and re-import it, but doesn't recognize the input. Has anyone experienced this problem or could provide a solution? Thank you!

> K562.hmm
GRanges object with 622257 ranges and 6 metadata columns:
seqnames ranges strand | V4 V5
|
1 chr1 10000-10600 * | 15_Repetitive/CNV 0
2 chr1 10600-10937 * | 13_Heterochrom/lo 0
3 chr1 10937-11937 * | 8_Insulator 0
4 chr1 11937-12337 * | 5_Strong_Enhancer 0
5 chr1 12337-13137 * | 7_Weak_Enhancer 0
... ... ... ... . ... ...
622253 chrX 155256806-155257806 * | 11_Weak_Txn 0
622254 chrX 155257806-155259206 * | 8_Insulator 0
622255 chrX 155259206-155259406 * | 6_Weak_Enhancer 0
622256 chrX 155259406-155259606 * | 7_Weak_Enhancer 0
622257 chrX 155259606-155260406 * | 15_Repetitive/CNV 0
V6 V7 V8 V9

1 . 10000 10600 245,245,245
2 . 10600 10937 245,245,245
3 . 10937 11937 10,190,254
4 . 11937 12337 250,202,0
5 . 12337 13137 255,252,4
... ... ... ... ...
622253 . 155256806 155257806 153,255,102
622254 . 155257806 155259206 10,190,254
622255 . 155259206 155259406 255,252,4
622256 . 155259406 155259606 255,252,4
622257 . 155259606 155260406 245,245,245
-------
seqinfo: 23 sequences from an unspecified genome; no seqlengths
> kpPlotRegions(kp, K562.hmm, col=K562.hmm$V9, r0=0.22, r1=0.3)
Error in (function (xleft, ybottom, xright, ytop, density = NULL, angle = 45, :
Unzulässige Farbspezifikation "245,245,245"

Answer: The rgb conversion shown [here](https://stackoverflow.com/questions/31574480/rgb-to-hex-converter) should help. I.e.

```
sapply(strsplit(K562.hmm$V9, ","), function(x)
rgb(x[1], x[2], x[3], maxColorValue=255))
```
Question: Hi all.

I am working on a matched tumor-normal somatic variant calling pipeline.
My pipeline is as below:

1. bwa mem alignment
2. sort and deduplicate with samtools and picard
3. Realignment with GATK
4. Base recalibration with GATK
5. Somatic mutation calling with Mutect / Varscan2
(For MuTect, I try both 1.1.7 and 2, with default setting. For VarScan2, I filter reads with mapping quality 20 and use processSomatic to pick out high-confidence calls)

Here are my questions:

1. When I filtered out the variants in introns/UTR/ncRNA, there are very little of intersection between Mutect/VarScan hit. The intersection between Mutect and Mutect2 is also very low. I am aware of the fact that the false positive rate is very high in current somatic mutation calling tools, but is there a way (a combination of parameter setting) that can filter out most of noises? (I know MuTect2 gives INDEL calling while old ones don't.)

2. I try to find a gold standard reference for whole exome sequencing. But what I've found so far are some articles using NA12878, simulating tumor mutation based on normal sample. Is there any reference I can use to evaluate my pipeline?

3. COLO829 is another candidate for me to use as reference.
Since it is a genome sequencing sample, would it be ideal to use it as a reference standard, by using the exonic intervals?

I am still a novice in WES. Any reply would be greatly appreciated.

Thanks.

Answer: The DREAM challenge consortium has generated [synthetic cancer data sets][1] for benchmarking (whole genome, but you could easily filter to WES after alignment). Brad Chapman et al. at Blue Collar Bioinformatics have validated a lot of mutation-calling tools against this data (see [here][2]).

[1]: https://www.synapse.org/#!Synapse:syn312572/wiki/62018
[2]: http://bcb.io/2015/03/05/cancerval/
Question: Hi everyone,

First post here - this forum is being of immense help in my bioinformatics journey so far.

Brief explanation: I have Illumina MiSeq 2x300 bp reads metagenomes. I have used BBduk + Trimmomatic to remove adapters and to quality trim the sequences. I have four output files - forward paired, forward unpaired, reverse paired and reverse unpaired. I run FastQC on all of them and the quality of unpaired output is slightly worse than that of the paired output.

Input Read Pairs: 3163058 Both Surviving: 2631476 (83.19%) Forward Only Surviving: 363260 (11.48%) Reverse Only Surviving: 48940 (1.55%) Dropped: 119382 (3.77%)

Forward only surviving: % of forward reads that was high quality but couldn't be kept because the paired reverse read was low quality
Reverse only surviving: % of reverse reads that was high quality but couldn't be kept because the paired forward read was low quality
Dropped: sequences which have been dropped because BOTH forward and reverse were bad quality.

From what I understood, Trimmomatic drops both forward and reverse reads when one or both of the reads do not go through the quality threshold. The whole pair will be dropped and will end in the unpaired output.

First question: why is that? Even though the quality is not high enough, the sequences will BOTH be dropped and end in the UNPAIRED output, but they will still be paired, right? So why is it called unpaired?

Second question: should I include my unpaired output in the assembly process? In my opinion, it would add a lot of information. I am going to use Megahit, and I would run it in a way that includes both paired and unpaired sequences. Like this:

megahit -1 1_R1_paired, 1_R1_unpaired -2 1_R2_paired, 1_R2_unpaired -o output/

Is it Megahit going to recognize the R1_unpaired and the R2_unpaired as still paired?

Sorry for any confusion I might have created, if necessary I will try to explain it better. Thanks in advance for your help.

Stefano

Answer: If I remember correctly: if you apply a lengthfilter-step and one of the mates is too short after the trimmers/clippers are applied, then this read will be dropped (i.e. it will never end up in either output) and its mate will end up in either unpaired forward or unpaired reverse.

Second, if your assembler allows it, you might as well add the unpaired reads. They might potentially cover some extra sequence or just bump coverage over the edge. In my experience, it doesn't matter whether you include them or not, but that is anecdotal and personal preference. What you lose is the paired-information, so the unpaired reads cannot help in the contigging/scaffolding process.

Edit: grammar/clarity
Question: Hello,

I will analyze soon RNA-Seq samples build using a TruSeq RNA Exome (illumina) from FFPE. I've already a standard RNA-Seq pipleline that works well on fresh tissue ( RiboZero kits). Pipeline is a classic STAR-DESeq2 + additional post-hoc analysis (splicing, fusion, pathways, etc..).

In this context, should I be carefull on some steps using capture-based RNA-Seq samples from FFPE samples ?

thank you

Answer: I would be interested to know if there are some criticals bias too !

A start in this Intro : https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5976456/
Question: After I got the filtered vcf file of snps with gatk pipeline, I tried to run PCA with gcta.
afterwords, I tried to find the explained variation percentage for each PCs(principal components).
what I'm confused is that the eigenvalue results I get are always the same number as the sample size.
I thought that the amount of possible PCs were the same as the variable dimension I have (which is the number of snps called, about 17000000 in this case) and there should be the same amount of eigenvalues related to it. Of course there could be same eigenvalues and PCs that explain small variation are useless, but isn't it possible that there could be 17000000 PCs and eigenvalues? So I thought when I want to get the explained variation percentage of PC1, I had to divide it by the some of 17000000 eigenvalues. could someone explain why this is wrong and why I always get the same number of eigenvalues as the sample size?

Answer: > I thought that the amount of possible PCs were the same as the variable dimension

If you think about it, if you have 2 points in a 3-dimensional space, you can rotate and shift the original axes in such a way that one axis passes exactly through the 2 datapoints. Therefore you don't need three principal components, with just 1 PC you can describe the 2 datapoints. In fact, if you have more dimensions than datapoints the number of PCs is at most the number of datapoints - 1. See also https://stats.stackexchange.com/questions/123318/why-are-there-only-n-1-principal-components-for-n-data-if-the-number-of-dime .

However, R's prcomp will return as many PCs as datapoints with the last PC with having stdev very close to 0. I'm not sure if this is because of the underlying linear algebra of SVD or because of the imprecision of floating-point numbers (*I think* the latter).
Question: Hello! I am trying to run RepeatMasker on some genomic sequences from some various teleost fishes. RepeatMasker starts and gives me this:

RepeatMasker version 4.1.1
Search Engine: NCBI/RMBLAST [ 2.10.0+ ]

Using Master RepeatMasker Database:
...path/bin/miniconda3/share/RepeatMasker/Libraries/RepeatMaskerLib.h5
Title :
Version :
Date :
Families :

Species "teleost fish" is not known to RepeatMasker. There may
not be any TE families defined in the libraries for this
species/clade or there may be an error in the spelling.
Please check your entry against the NCBI Taxonomy database
and/or try using a broader clade or related species instead.
The full list of species/clades defined in the library may be
obtained using the famdb.py script.

It does the same with any other species I've tried including danio and homo sapiens so I think the issue is the database. Please help! Thanks.

Answer: I figured it out!!! RepeatMasker is running python3, so when I was running just python it was finding h5py but when I ran python3 it couldn't find it.

I fixed it by adding:

sys.path.append("/mypath/bin/miniconda3/lib/python3.7/site-packages")

directly to the famdb.py file, between import sys and import h5py

and it worked!! Hope this helps someone else