cannin/biostars_qa huggingface dataset:
Question: Hello everyone,

I want to explore differential expression of a relatively small set of genes (~30). I used this set as reference to map reads.

I used Bowtie2 to build an index of these genes sequences, then I mapped transcriptomes reads to this reference, and I got the sam files. In the sam file, for each mapped reads, there is the gene name to which the read mapped as RNAME (reference sequence name).

Now, I want to obtain the matrix to see if there are differentially expressed genes.
Looking to HTSeq, I noticed he wants gff or gtf file to build this matrix, but since I don't have the whole genome as reference, I cannot use it.

Do you have some suggestion on how to build the matrix? Should I use samtools and build a bed file to get the coverage for each gene?

Answer: There are a couple of problems with this appoach:

1. `bowtie2` is not a splice-aware aligner so reads spanning exon/intron junctions will likely go unmapped.
2. Alignment against a subset of the total genome or transcriptome is problematic as the aligner will always try to find best matches for every read. If the true origin is a gene that is not included, you will get false-positives (maybe even a lot of them).
3. Most differential analysis frameworks rely on normalization based on a large number of genes not changing between conditions, e.g. TMM from `edgeR`. With only 30 genes you will have a hard time doing proper normalization. Don't use RPKM or FPKM as they do not correct for library composition changes, e.g. where a few genes with high counts in one condition "dominate" the sequencing run (=take away reads from other genes leading to fewer counts even though these are not differentially expressed).
4. How do you identify differentially-expressed genes (DEG)? I am pretty sure established tools do not work on such small subsets. I hope you do not plan custom statistics such as pairwise rank sum tests or t-tests?

I suggest you align against the entire genome with a splice-aware aligner such as `hisat2` or `star` or use a lightweight quantifier such as `salmon` or `kallisto`. Then use established DEG tools such as `DESeq2`, `edgeR` or `limma` to get all DE genes and check how many of ones you are interested in are among them.
Question: I have few VCF Files which contains Chromosome number, position, Variation allele (Ref and Alt).

How can I generate the Gene name, its Refseq and listing all of its transcripts?
1.	Any tool which gives these kind of results when upload VCF data?
2.	Is there any module in Biopython or Bioconductor or script for this?
3.	Is there any lists that links gene name to Refseq and transcript information?

Thank you for the valuable suggestions

Answer: Ensembl BioMart is probably your best bet.

You may need to use/write a quick script or awk command to extract your chromosome number and position but you can enter a list of chromosomal positions in the filter menu under regions->multiple chromosomal positions. You can then choose which attributes you would like to download such as Ensembl Gene ID or (under the external references menu) Refseq id.

If your list is too long you'll probably have to use the Perl API. Short tutorial available here: http://www.ensembl.org/info/data/biomart/biomart_perl_api.html
Question: Hi all,
Does anyone know how to make database for compressed file fasta.gz using blast? I am using the blast2.2.29+. I saw the tutorial of makeblastdb, there is no option for dealing with the fasta.gz file. Does that mean I have to unzip the file first and then build database?

Answer: You should be able to uncompress the file then pipe the output to makeblastdb, without an intermediate file. The command would look something like:

gunzip -c myfastafile.gz | makeblastdb -in stdin <other options>

or

gunzip -c myfastafile.gz | makeblastdb -in - <other options>
Question: Hi all, I am new to the BioInformatics, and quite a beginner in programming languages. Can anyone suggest me some sources where I can at least learn 50% of the scRNA seq data analysis? I am familiar with C language, and I know a little bit of molecular biology too.

Answer: Comprehensive collection of all things single cell (including tutorials) : https://github.com/seandavi/awesome-single-cell
Question: Hello dears all,

Actually I have a .frq file (allele frequency file) which It is generated from .vcf file through VCFtools and it is included by five populations (The vcf file is generated by Stacks pipeline). But the point is, I don't know how can I visual the allele frequency. I will be so grateful if anyone can help me to solve this issue. Here is a part of my .frq file.

Best regards.

CHROM	POS	N_ALLELES	N_CHR	{ALLELE:FREQ}
CM003279.1	116397	2	4	A:0.5	C:0.5
CM003279.1	184945	2	4	C:0.5	G:0.5
CM003279.1	185823	2	36	C:0.944444	T:0.0555556
CM003279.1	408595	2	54	T:0.5	A:0.5
CM003279.1	420946	2	10	A:0.5	T:0.5
CM003279.1	420946	2	42	T:0.619048	A:0.380952
CM003279.1	432468	2	4	A:0.5	T:0.5
CM003279.1	432468	2	4	T:0.5	A:0.5
CM003279.1	521247	2	6	A:0.5	C:0.5
CM003279.1	521247	2	4	C:0.5	A:0.5
CM003279.1	768668	2	46	T:0.978261	A:0.0217391
CM003279.1	768678	2	50	A:0.84	G:0.16
CM003279.1	768679	2	50	T:0.64	C:0.36
CM003279.1	768685	2	50	G:0.84	C:0.16
CM003279.1	768691	2	50	G:0.84	A:0.16
CM003279.1	768705	2	50	C:0.84	T:0.16
CM003279.1	768709	2	54	A:0.814815	T:0.185185

Answer: Perhaps you can use online VEP (Variant Effector Prediction) server from ENSEMBL, just google it. I don't know whether your VCF file derived from a plant species or animal species or prokaryotic species like bacteria. Depends on the species, you can choose the correct VEP-ENSEMBL server. Here, you could predict variants effect (with SIFT score) for each position of the SNPs from your VCF input file and you have to choose the correct species of interest when you give your input VCF file. For instance, here is the VEP ENSEMBL server for plants: https://plants.ensembl.org/Triticum_aestivum/Tools/VEP?db=core you can read here for variant effects prediction/variant of alleles and annotations concerning VEP server: https://plants.ensembl.org/Help/View?id=484 Hope it may help.
Question: Hello,
we're getting RNAseq data from increasing numbers of "emerging" model species for which little is known. The data often wasn't meant to provide population genetic estimates, but could perhaps be used to provide some.

Is it possible to estimate Ne from the SNP patterns that can be identified from:

a single. sequence pooled from many individuals from a single population?

a single diploid genome?

RNAseq data where multiple individual from a single population (siblings or not) were independently sequenced.

Cheers,
yannick

A somewhat related question, but specific to pools is here.

Answer: It seems that people use RNA-seq or exome capture data for demographic inference by using synonymous mutations, thereby getting around the issue of non-synonymous mutations being non-neutral. Therefore, it should be similarly possible to use your RNA-seq data to estimate Ne (your third point above) if you use the synonymous mutations only?

For example, see:
> Fraïsse C, Roux C, Gagnaire P, Romiguier J, Faivre N, Welch JJ, Bierne N. (2018) The divergence history of European blue mussel species reconstructed from Approximate Bayesian Computation: the effects of sequencing techniques and sampling strategies. PeerJ 6:e5198 https://doi.org/10.7717/peerj.5198
Question: I have 1000 fasta files that have simulated reads, and I want to split each of these 1000 files into separate files (one per chromosome) as I need this for some further analysis. The header lines look like this >chr1:startpos:endpos.

I wrote this code (https://gist.github.com/ethanagbaker/6e40c58127b7ca8b9242) in Python, and it works, but it is very slow (ie it has been running for more than 24 hours on a cluster). This seems like it should be a really quick thing to do. Is there a better way to be doing this?

Thanks!

Answer: $ pip install pyfaidx
$ faidx -x sequences.fa

Or you can use the `Fasta` class and write your own script to do the same thing. Your code is slow because it is opening a bunch of files in a loop, and then opening (the same files?) and reading them completely to get the sequences. Indexed file access will be orders of magnitude faster.
Question: let's say I have a fasta of a protein sequence

> albumin

MKWVTFISLL FLFSSAYSRG ... ... ...

I want to split the sequence into all possible consecutive 8 amino acids (only in 1 direction, amino -> carboxyl) (And no looping(I don't know if it is the right expression), say, GMKWVTFIS)

I need

> fasta.albumin1

MKWVTFIS

>fasta.albumin2

KWVTFISL

> fasta.albumin3

WVTFISLLF

...

> fasta.albumin13

FSSAYSRG

And, I want to do this for all known human protein sequences.

How would I do it???
I need the result as a fasta or fasta files. And IDs for resulting 8-mer seuqeunces need to be unique.

Answer: Using [bioalcidaejdk][1]:

$ java -jar dist/bioalcidaejdk.jar -e 'final int N=8;stream().filter(F->F.length()>=N).forEach(F->IntStream.range(0,F.length()-N).forEach(L->{out.println(">"+F.getName()+"."+(L+1)+"\n"+F.subSequence(L,L+N));}));' input.fasta

[1]: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html
Question: I know that `blastp` has a `-remote` option, which allows you to run your query on NCBI servers (see [here](https://www.ncbi.nlm.nih.gov/books/NBK279680/)).

Is there a way I can set up my own remote server that could be queried using `blastp` and would run queries against my own custom database? I am running lots of `blastp` queries on a cluster of machines, and I think it would be much faster if I delegated this work to a single machine that could use a local blast database, instead of querying a database stored on a distributed file system?

Answer: Check the [documentation for NCBI BLAST AMI][1]. While this is talking about enabling remote blast on AMI, the idea should be similar for local implementation.

[1]: https://ncbi.github.io/blast-cloud/doc/setting-up-blast.html
Question: I want to uninstall *Bioconductor* [affycoretools package](https://www.bioconductor.org/packages/release/bioc/html/affycoretools.html) and reinstall it, who can help me?

Answer: Go to the **Packages** in right bottom corner of Rstudio, sear the package name and click on the adjacent **X** icon to remove it.

![Rstudio Screenshot][1]

[1]: https://image.ibb.co/j4KZjG/test.png

Reinstall the package from Bioconductor/CRAN.

##Bioconductor
source("https://bioconductor.org/biocLite.R")
biocLite("affycoretools")
Question: Hello everyone,

Today I spent 4 hours on a python script for executing a shell command constructed using user inputs within the script, still It was unsatisfactory to what I wanted. I achieved same thing using R < 30 mins. I am in under obsession that one of the scripting language is very important. I use R, shell and terminal utilities for most of my tasks.

So, the argument would be I haven't come across to such a task in exploratory next generation sequencing data analysis, where you specifically need a scripting language like Perl/Python or Ruby.

Thus, if one is not writing new tools from scratch, how important is to learn any of the scripting language. I used Perl in start and then got influenced by python co-workers that python is better and more extensive and now lost in both. But I am confident in R.

Could someone also elaborate, the part and parcel of work, where there is no move without any of the scripting language. One can also argue about the audience, perl user , python and R users at the consumer end, if you make software but what if the results and tasks are for oneself or its just a matter of taste.

Thanks for your input

Answer: From a certain point of view you are right: you don't need to learn python or perl programming, and you can do a lot of work just in R. CRAN and Bioconducator contain a lot of libraries, and R allows you to create good plots, which is one of the most important skills for a bioinformatician. I know people who mostly use only R, and they are very good in their work, so you don't have to be afraid if you don't want to learn other programming languages.

However, consider that learning at least the basics of other languages does not take too much time, adds skills to your curriculum, and allows you to learn new approaches to programming. R is a language for data analysis, centered on the data.frame structure, and is suited mostly to analyze data organized in tabular form. That's a good paradigm for analysing data, but as soon as you have to other problems R becomes more clunky and less useful. Web programming is an example, but also data manipulation (convert one format to another, automatizing tasks on the shell, etc.) is important. If you know the philosophy under which other programming languages are developed, and what tasks are they typically used for, your programming skills will be stronger, and you will be able to choose among a larger range of approaches when you have to solve a certain programming topic.
Question: Hello!

I'm sorry if this question comes of as naive or ignorant because I'm very new to Bioinformatics. I'm trying to do an alignment with STAR and was wondering if I could access a pre-made STAR index for the mm10 genome. I was told I could do this from UCSC but have had no luck finding it there.

So my question is, Are there pre-made STAR index files for the mm10 genome that I could download? And if so where and how?

Thanks in advance for any help and I'm sorry to ask such a trivial question! Let me know if there's anymore detail I can give!

Answer: You can generate indexes yourself easily enough. Follow the directions here: https://www.biostars.org/p/198315/ . MM10 genome [from UCSC is here][1].

@Alex has some pre-made indexes available at [STAR Genomes site.][2] There does not appear to be a UCSC version of Mouse but there is [Gencode Mouse][3] which you can use.

[1]: http://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz
[2]: http://labshare.cshl.edu/shares/gingeraslab/www-data/dobin/STAR/STARgenomes/
[3]: http://labshare.cshl.edu/shares/gingeraslab/www-data/dobin/STAR/STARgenomes/GENCODE/GRCm38_GencodeM8/
Question: Hello,

I am running GATK 3.8 on a Linux cluster. I am trying to run the tool VariantFiltration with the following commands:

java -jar $GATK/GenomeAnalysisTK.jar \
-T VariantFiltration \
-R $ref \
-V $wd/vcf/filtered/stick84_SNPs_GATK.vcf \
--filterExpression "QD < 2.0" --filterName QD2 \
--filterExpression "FS > 60.0" --filterName FS60 \
--filterExpression "MQ < 40.0" --filterName MQ40 \
--filterExpression "MQRankSum < -12.5" --filterName MQRS-12.5 \
--filterExpression "ReadPosRankSum < -8.0" --filterName RPRS-8 \
-o $wd/vcf/filtered/stick84_SNPs_filtered_GATK.vcf

However, it is throwing me this error:

ERROR MESSAGE: For input string: "nan"

I don't know what this error means and how to fix it. All the potential help I find online sends me to ghost pages of GATK that I cannot access :(

Please, any help will be greatly appreciated! Thank you and take care!

Answer: Ah, I've found out the reason for that here: https://github.com/broadinstitute/gatk/issues/5582. I hope this helps others as lost as me! I did not realize the INFO fields could have "NaN" values. Solution to fix is:

bcftools view in.vcf.gz |
sed 's/=nan/=NaN/g' |
bgzip > fixed.vcf.gz
Question: There exists a lot of literature on distinguishing driver mutations from passengers. I am trying to build my own deep learning model to do the same. I am facing some potential issues. First, I have downloaded COSMIC mutation data and used the FATHMM labels to designate drivers (positive examples) in my dataset. I am sceptical to use passengers from COSMIC as they may be false negatives. So I turned to the 1000 genome project to download SNVs (to construct my negative examples). I am unsure if this is correct, however, I have seen some papers do the same. Do I need to apply any filters on the 1000 genome SNV data to construct the final dataset? One such paper talks of using SNVs with a global minor allele frequency≤1%.

Answer: You'll fist need to clarify what exactly you mean by driver mutation vs passenger mutation. 1) Are you concerned with somatic mutations or germline mutations (most often it is somatic, if you say "driver mutation")? 2) Are you concerned mostly about protein-coding mutations or non-coding mutations?

1) If you are interested in somatic driver mutations, I highly recommend that you do NOT use germline variants from sources such as the 1000 genome project as "passengers". Germline variants are systematically different than somatic mutations because they undergo substantially more negative selection than somatic mutations do (PMID: 19654296). If you are interested in "pathogenic" germline variants that relate to cancer, then using high allele frequency germline variants as "passengers" is an ok but not perfect solution (see Kevin's comments).

2) If you are interested in protein-coding driver mutations, I highly recommend that you do not develop a method that also tries to predict non-coding driver mutations. The recent PCAWG papers have highlighted that indeed most of the somatic driver mutations in cancer happen in protein-coding regions. Moreover, non-coding prediction methods applied to protein coding regions generally don't fare well in benchmarks (see PMID: 32079540). The top methods in that recent benchmark at predicting driver mutations were CHASM, CTAT-cancer (see my paper, PMID: 29625053), DEOGEN2 and PrimateAI. I recently released CHASMplus (pmid: 31202631) which is a substantial improvement over top performing CHASM at predicting driver mutations, and my training labels are available (https://chasmplus.readthedocs.io/en/latest/faq.html , see "Where can I obtain the training data for CHASMplus?").
Question: Hi All,

can anyone recommend me a tool/package, which allows me to create Venn diagrams up to 11 sets? The packages which I have found so far can support to create less sets only.

Many thanks!

Answer: This is not an answer, and you may be upset :-) In my opinion, it is not practically useful to do Venn diagram for more then 4 or 5 sets. More than that, it is horribly difficult for reader to look at and understand what it conveys. It is also difficult to explain about the result for the authors themselves. So, if possible, divide them into smaller groups. I beg you as a reader.
Question: I have a long list of complete bacterial organism names (more than 100000, thus impossible to search and download it line by line). Format is one name on each line. I need to download GCA (It must be GCA instead of GCF) fasta files of the corresponding genomes from https://www.ncbi.nlm.nih.gov/genome/browse/ (Specify Levels as Complete).

I have to achieve this through command lines. How to do it efficiently? Thank you.

Answer: cat species.txt
Porphyromonas levii
Porphyromonas somerae

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/bacteria/assembly_summary.txt

IFS=$'\n'; for next in $(cat species.txt); do awk -v SPECIES=^"$next" 'BEGIN{FS="\t"}{if($8 ~ SPECIES && $12=="Complete Genome"){print $20}}' assembly_summary.txt \
| awk 'BEGIN{OFS=FS="/"}{print "wget "$0,$NF"_genomic.fna.gz"}'; done \
| sh

NOTE: Only 8,413 Bacterial genomes have "Complete Genome" assembly level status (not even 10% of your list of names). For example, nothing will be downloaded in the example shown above. Do you really need to limit yourself to such a small subset?

1577 Chromosome
8413 Complete Genome
52594 Contig
54565 Scaffold
Question: I need to do some statistical calculation on CNV data that is publicly available in TCGA website. I am from computer background and no idea about the meaning of these data. I tried searching for the meaning of the files, how they are structured and all but nothing was helpful. It would be nice if someone could give a short overview of the meaning of data in each column. A sample of data present in a CNV file

Sample	Chromosome	Start	End	Num_Probes	Segment_Mean
DEBUT_p_TCGAb45_81_wRedosSNP_N_GenomeWideSNP_6_G03_729530	1	61735	415164	28	-0.0504
DEBUT_p_TCGAb45_81_wRedosSNP_N_GenomeWideSNP_6_G03_729530	1	462793	629241	4	1.822
DEBUT_p_TCGAb45_81_wRedosSNP_N_GenomeWideSNP_6_G03_729530	1	668210	2138242	350	-0.0311

So I want to know what is the underlying meaning of sample, how they are named. Sample name seems always a big name. Then I understand the columns chromosome, start and end. But I don't understand the meaning of Num_probes and Segmentation_mean?

Answer: Sample should be the unique sample identifier for this SNP-array experiment, it corresponds to a unique sample run by TCGA and can be a normal or tumor sample. Will map to a unique id like this: https://wiki.nci.nih.gov/display/TCGA/TCGA+barcode

The rows correspond to contiguous chunks along the chromosome with the same DNA copy-number. "Chromosome" is the chromosome, can be 1-22, X or Y (see human genome). Start is the physical start location for the segment along said linear chromosome, end is the end coordinate. Num_probes is the number of SNP-array probes falling within the segment (these were used to call copy numbers). Segment_Mean is the estimated copy number for that particular segment.

So this sample has (log2) copynumber -0.0504 on chromosome 1 from bases 61735 to 417164, at 462793 along the same chromosome the copy number changes to 1.822 (i.e. a genomic gain) and this segment ends at 629241. The next segment for which you have information is 668210 to 2138242 where the copy number has returned to a basically normal level. These segments should keep coming for chr1 until about 255M after which the chr2-segments should start.

The data should basically be log2( intensity of sample / intensity of reference with normal copynumber ), this means that segments with a normal diploid copy number are around log2( 2/2 )=0, single copy losses are at log2(1/2)=-1, homozygous deletions at log2(0/2)=-Inf, correspondingly for gains you go from log2(3/2) and upward with experimental noise leading to deviations from these discrete levels.

For some light reading about the concept of copy-number analysis check out these old papers: https://www.ncbi.nlm.nih.gov/pubmed/16899659 and https://genomebiology.biomedcentral.com/articles/10.1186/gb-2008-9-9-r136 and maybe http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-8-382
Question: Hello everyone!

I am trying to install MUMmer 3.9.4 on Ubuntu, following the instructions given here: https://github.com/gmarcais/mummer/blob/master/INSTALL.md

When I try to run nucmer, I get this error message:

> nucmer: error while loading shared libraries: libumdmummer.so.0:
> cannot open shared object file: No such file or directory

I have tried to look it up online, but can't seem to find anything.

Any help would be greatly appreciated!

Nancy

Answer: @nancydong20
I had the same problem as you, I did these steps and it helped:

1. This is critical: Uninstall the existing MUMmer with Synaptic Package Manager. (If you install with `sudo apt-get install mummer`you get v 3.22 and it is still called by any software which you use. This is why you cannot use the option *CPPFLAGS="-O3 -DSIXTYFOURBITS"* to install 64-bit Mummer).
2. Download the stable MUMer distro 3.23 from https://github.com/gmarcais/mummer
3. Extract the package
4. Install MUMer with this code:
5. `cd MUMmer3.23`
6. `mkdir aux_bin`
7. `make check`
8. `make install CPPFLAGS="-O3 -DSIXTYFOURBITS"`
9. ` mummerpath=`pwd` `
10. `cd ..`
11. `export PATH=$mergerpath:$mummerpath:$PATH`

In the end you will get:

- (maximum reference length is 2305843009213693948)
- (maximum query length is 18446744073709551615)
- process 6420479 characters per dot

My OS Ubuntu 14, 64 bit.
Hope it helps.