cannin/biostars_qa huggingface dataset:
Question: Hello everyone!
I'm new to bioinformatics, and I'm having a really hard time trying to make this work.
What I'm trying to set up is this https://github.com/NYU-BFX/lncRNA-screen

So I'm working with Long non-coding RNAs, and this pipeline created by Applied Bioinformatics Laboratories (New York, NY), does exactly what I need.
However, I'm finding quite hard to set it up, could anyone help me?

*It says it uses SGE which I only got it to work with docker, is SGE really necessary? I only have 1 machine.

*Needs to install and set to path r/3.3.0, python/2.7.3, java/1.8 and samtools/1.3

*It has a linked folder for my RNA-seq and Chip-seq but I don't know how that works.

*Also says I need https://github.com/NYU-BFX/RNA-Seq_Standard even if I have my own RNA-seq (which I do have).

*The documentation says sratoolkit is included, but, my lack of experience makes me not understand how that works.
*Here's a requirement file https://github.com/NYU-BFX/lncRNA-screen/blob/master/inputs/system_requirement.txt

This is my first post here, so I may do something wrong or post this question in the wrong place.

Answer: SGE is a scheduling/job submission system for computing clusters. You don't need it to run locally, though your machine better be a beast, as STAR uses a lot of RAM and is slow without several processors (as is every aligner). If you have only a few samples, you can probably get away with it, but if you have dozens, you're going to be waiting a while. I'd see if your organization has a computing cluster that you can get access to.

As for R, python, java, and samtools, they are all easy to install and add to your PATH. You can google how to do it depending on your system, and many distros can install them through package managers. Or you can look into Anaconda, which makes install all of those and automatically adding them to PATH very easy regardless of your OS. sratoolkit is also simple to install.

The links are symlinks, basically saying where your folders containing your ChIP-seq and RNA-seq data should be relative to that folder. In this case, it looks like you should have folders for them a level up from the installation directory. You can also just replace the link so that it points to wherever your data files are for each.

Installing something like this is a headache even for experienced bioinformaticists - relative lack of documentation, heavy reliance on relative paths, etc. I imagine it is one of those things that will result in a million errors with uninformative tracebacks that you'll spend days fixing before getting it to run in full. If you don't know how to install/add basic programs to your PATH, I would take a few days to learn how to do that and utilize the command line properly. Otherwise, you will likely continue to be frustrated.
Question: Hello,
Does anyone knows how to download or get only snp from hapmap EUR population for gtex v8 models?
I tried to get the snp annotation file from gtex portal and downloaded the lookup table in order to get the snp_annotation file.
Thank you.

Answer: Hi Rheab,

Yes, certainly! First, for JUST hapmap, go to [hapmap's FTP site][1].

One thing, though. Nowadays, there are more comprehensive options out there. Two that I recommend, for various reasons, are:

1. Check out [Ensembl VEP][2]. Steps: click link, paste in an rsID > run job > once job is finished, you will see a table in your browser. You can download this, including columns broken out by global MAF, then MAF for each ancestry specifically. [Example result for a variant I just ran][3].

2. If you want to take a deeper dive into a confined genomic locus (max 50,000bp in length) then [Kaviar][4] is a good option. This database contains more variants than dbSNP, even, and includes allele fractions down to very low percentages.

[1]: https://www.sanger.ac.uk/resources/downloads/human/hapmap3.html
[2]: https://www.ensembl.org/Tools/VEP
[3]: https://www.ensembl.org/Homo_sapiens/Tools/VEP/Results?tl=TykABz8JFAblrPFQ-8512135
[4]: https://db.systemsbiology.net/kaviar/
Question: I am working on the Differential Expression Analysis on RNA-Seq Data between tumor sample and normal sample, while my lecturer suggest me better use the raw expression data instead of using FPKM as expression unit , anyone could tell me why ? or what is the advantage for this ?

Answer: The short answer is...the statisticians who wrote the software to do differential expression analysis have thought a lot about how best to normalize, and you should just let them do it. Those softwares all want raw counts as input so they can normalize, you should give them that, and let them do their thing. If you know what you are doing, and have compelling reasons to apply a different normalization function, the softwares have ways to allow you to do that, but a simple FPKM method is not going to be better than what the authors of the software have worked out.
Question: I'm just curious if anybody can suggest command line programs that are designed to be truly parallelized?

I found a great post on running serial programs in parallel with GNU parallel: Tool: Gnu Parallel - Parallelize Serial Command Line Programs Without Changing Them

However, instead of using tricks to run serial programs in a parallel fashion, is there anything out there that can be run in parallel directly?

I'm mainly interested in programs used for chip-seq and rna-seq workflows. For example, I think STAR aligner can handle parallel processes, but I believe common tools like tophat, MACS, and cufflinks are all serial at some point.

Thanks in advance for any input.

Answer: Since you asked for tools for the analysis of RNA-seq and ChIP-seq data: there's [deepTools][1], [seqminer][2] and [spark][3] that all allow you to make use of multiple processors. I only have worked with deepTools though and of course, I think they're pretty great and useful (I helped develop them)

[1]: http://deeptools.github.io/
[2]: http://bips.u-strasbg.fr/seqminer/
[3]: http://www.sparkinsight.org/
Question: Hi,

I'd like to calculate the increase in qualitative disease score per risk allele. I've a cohort of 1000, and a single SNP, for each of the individuals, I've a score from 0 to 100. It was recommended that I use a linear regression model.

Can anyone elaborate on why you might use a linear regression model for this? Are there any other models that would fit this task? Suggestions on R packages I could use, are welcome.

Thanks

Answer: The versions of linear regression (well, a linear model really) that you're familiar with are ANOVAs and T-tests. This is essentially an ANOVA, possibly with a partial factorial design (depending on how you want to set things up).

BTW, the alternative to this would be to divide the scores by 100 and use logistic regression. This would be to get around some ceiling and floor effects in more complex models. Of course, if you only have essentially two alleles then a Wilcoxon test (aka Mann-Whitney U-test) could work. But that only makes sense if you just have two groups. That's essentially just a non-parametric T-test.

Anyway, graph the data and have a look. Then fit it with lm() in R and plot that too to see if the fit seems reasonable.
Question: Hello, everyone !

I'm looking for a way to use HTseq with multi-thread.

I couldn't find any options about multi-thread.

Anybody knows how to, please ?

(I know there are tools support multi-thread like STAR, HISAT2. but just wonder whether HTseq doesn't support it.)

Thank you. : )

Answer: HTSeq-count (I guess that is what you are referring to) doesn`t support it as of now for a single file, see the issue for the discussion and alternatives: https://github.com/htseq/htseq/issues/7#issuecomment-638075405

To be precise (citing the developer)
> That's it. The -n can process up to n BAM files in parallel to speed up, but cannot use multiple cores on the same file.

Other HTSeq commands are muti-threaded and support the `-n` option though.

You could use featureCounts instead if you really need multi-threading support. Or run HTSeq count through `parallel` on multiple files, but this will not scale favorably complexity-wise (duplicated coordinate parsing, duplicated tree structures, increased simultaneous file IO), so this is not recommended.
Question: Hi,
I was taking a lot at the variant allele frequencies in TCGA NSCLC data (from the the MAF somatic SNV file). I have seen thet the highest VAF is 50% and it apears that they have a cutoff there. Does it mean that a SNV of 60% vaf can not be a homozygous mutations in 60% of cells (while tumor purity of 60%) or a subclonal snv at higher purities?

As well there is a lot of variants with 0% VAFs! which is also strange.

What I want to learn from this story is: what should I consider when facing a oncoplot? can a variant on oncoplot have only 3% vaf? specially in paired primary vs recurrent tumors, can it be that a variant is simply missing in one because it has a less than cutoff level vaf ?

Best

Answer: Answering myself:
because mutations are considered to happen in heterozygous state, therefore at 100% purity you will have a 50% VAF. VAF of more than 50% belongs to the LOHs.
Question: dsk readme shows command line with -file and -kmer-size.

The one I used (no --version, so don't know which) on genouest.org cluster doesn't understand those options.
It seems to be only dsk file_in kmer_size.

So time for you to update the readme or for me to update dsk ?

Answer: Sebastien,

You need to update your DSK version.
[https://github.com/GATB/dsk/releases][1]

[1]: https://github.com/GATB/dsk/releases
Question: Hi,

I have a vcf file which I would like to filter by a list of SNP positions, which are present in a positions.txt file.

bcftools filter -i 'POS=@positions.txt' input.vcf.bgz > output.vcf

However, this says:

Error: could not parse the expression. Note that the "@file_name" syntax can be currently used with ID column only.

Is there another simple way filter by POS in bcftools?

Thanks.

Answer: Yeah, I figured.

Actually found the answer here, so you can close the thread

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

thanks.
Question: I am using STAR aligner to align an RNA-seq data with a human reference genome and I was analysing the log.final.out file. I am not sure what the following results mean:

Number of splices: GT/AG

Number of splices: GC/AG

Number of splices: AT/AC

Number of splices: Non-canonical :

I understand that the number of splices (total) gives the total number of splices in all the reads and number of splices (annotated) gives the those which are present in the annotation file.

What does canonical/non-canonical splicing mean?

Answer: The canonical splice sites are those originally described and most commonly found (like in ~99% of introns) and have GT at the donor site (just after the 5' end of the cut) and AG at the acceptor site (just before the 3' end of the cut). The rarer variants GC/AG and AT/AC are also considered canonical and all other dinucleotide combinations are non-canonical.
Question: I have a bash script about HaplotypeCaller.sh. Is it possible to run this in gnu parallel ?

cat HaplotypeCaller.sh | parallel

#!/bin/bash

for i in *.bam; do
gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I $i -O $i.vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list; done

Answer: The answer given by [Malcolm.Cook][1], [ole.tange][2] and [cpad0112][3] works well

Method 1:

parallel gatk --java-options "-Xmx40g" HaplotypeCaller \
-R /media/gatk/Homo_sapiens_assembly38.fasta -I {} -O {.}.vcf.gz \
--dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz \
-L /media/gatk/target/ag_V6_list.interval_list ::: *.bam

Method 2:

#!/bin/bash

do_one() {
gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I "$1" -O "$1".vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list
}
export -f do_one
# test that do_one works on a single bam file. Then:
parallel do_one ::: *.bam

[1]: https://www.biostars.org/u/266/
[2]: https://www.biostars.org/u/6149/
[3]: https://www.biostars.org/u/20673/
Question: Is there a way to search all human or mouse protein sequences for a short (~5 amino acids) degenerate amino acid sequence?

Answer: You can do this easily by installing the FAST: Fast Analysis of Sequences Toolbox ([publication][1])([github][2]).

You can install it by using this command (only use sudo if you need to):

(sudo) perl -MCPAN -e 'install FAST'

You will also need the fasta files containing human and mouse proteins. You can get these from NCBI.

Lets say that the sequences you are looking for are `ASNNF`, `ASNLF`, or `ASNKF`. To search for this you can use fasgrep which is a utility of FAST:

fasgrep -is "ASN[NLK]F"

This will return sequences that have either `ASNNF`, `ASNLF`, or `ASNKF`.

[1]: http://journal.frontiersin.org/article/10.3389/fgene.2015.00172/abstract
[2]: https://github.com/tlawrence3/FAST
Question: I need to **pairwise align protein structures**.

Specifically, my task requires

- alignment of structure predictions in PDB format, for EACH of ~1000
full length proteins

**versus**

- X-ray crystallography based structure PDB files for Pfam/hmmsearch
delimited domain only regions of EACH of these 6 (shorter domain)
sequences.

This means **~1000*6 = 6,000 pairwise structure alignment runs**.

Which software do you recommend from [Wiki][1] or elsewhere, with my following requirements:

1. I should be able to download and run on my local laptop or
university HPC, and not on a webserver (unless batch submission and
download are allowed).
2. The runtime, diskspace and RAM needed for each of my 6,000 pairwise
alignment need to be reasonable for either my local laptop or on my univ's
HPC (with ~2GB RAM nodes and 12cpu limit on my user account)
3. I need to be able to parse the alignment results programmatically
(using Perl or Python scripting...)

**FINAL GOAL** = Classify and/or rank each of the full-length proteins based on likelihood of domain of interest being "present", **not** using pairwise **sequence alignments**, **but** from reported metrics of pairwise **structure alignments** (rather than my arbitrary yes/no/maybe classification)

If your suggested bioinformatics pipeline has been published, better yet!

Thank you all, in advance. Cheers!

[1]: https://en.wikipedia.org/wiki/Structural_alignment_software

Answer: These links will hopefully be useful as references:

- https://en.wikipedia.org/wiki/Structural_alignment_software
- https://en.wikipedia.org/wiki/Structural_alignment
- https://proteopedia.org/wiki/index.php/Structure_superposition_tools

I have tried locally TM-align, Fr-TM-align, Matt, FATCAT, MAMMOTH, SSAP, FAST, THESEUS, and probably couple of others that I have forgotten. All of them used to be available as independent binaries, but for some of them that may not be the case due to age or general lack of support. All of them are reasonably fast, so I don't think that would be a limiting factor for you no matter what you choose. They have various advantages and disadvantages, and also their outputs are different in terms of information provided. I suggest you try several programs and see what functionality is most useful for your purpose.
Question: I'm using Rvtests for association of exome chip genotypes with a quantitative phenotype. The SKAT test runs fine with default parameters, but when I try to change the alpha error and variant weighting parameters like this:

--kernel skat[nPerm=100:alpha=0.001:beta1=1:beta2=20]

I get an error message:

zsh: no matches found: skat[nPerm=100:alpha=0.001:beta1=1:beta2=20]

Any help would be greatly appreciated...

Thanks!

Luca

Answer: Problem solved... The zsh shell had trouble sorting out the `[]`. Switching to bash shell works.
Question: Hi Biostars,

I got two output files (*genes.results and *isoforms.results) from aligning reads to the Trinity denovo assembly(Trinity.fasta) with Bowtie2 and expected read counts with RSEM. Out of two output files, which one should I feed to DESeq2 for identifying DEGs.

Any suggestions will be appreciated. Thanks in advance.

Answer: Depends on what do you really want, the genes table represents the global expression of the gene (as a summary of all isoforms for the gene), the isoforms table is for each transcript of the gene. So, if you run DESeq2 in genes, you get genes, and isoforms, you get transcripts.

From a biological point of view, genes are representing a general function, the isoforms a more controlled regulation, if the biological response is to shut down the gene, you can expect all isoforms will be downregulated, but if the biological response requires to shut down a particular isoform, the gene will be not DE, but one isoform will be.
Question: Hi everyone,

I'm wondering what would be the best strategy to obtain the list of all k-mers in the graph with each k-mer labeled as reference (i.e., from specific path) or not-reference. A very nice to have extra feature would be to have, for each k-mer, the index(es) of its occurrence(s) in the reference (linear) space.

I see two possible roads:

vg kmers ...
vg find -k ...

Both strategies have prons and cons. The former produces more succinct output but I haven't found a way to discriminate whether k-mers belongs or not to a path (nor I found an easy way to obtain the index). The latter has more extensive output, but also produces much larger outputs.

Thanks for any suggestion,
Michele S.

Answer: The easiest way is probably doing it outside vg:

1. Extract all kmers (kmer occurrences) with `vg kmers`.
2. Extract the selected paths in FASTA format with `vg paths -v graph.vg -F -p path-names.txt > output.fa`.
3. Determine the kmers in the paths using an external tool.
4. Compare the kmer sets with an external tool.

The path positions for non-path kmers are not always well-defined. There is some machinery for determining them (for alignments, not kmers), but even then, we often have to make arbitrary choices or give up trying.
Question: I have a bash script about HaplotypeCaller.sh. Is it possible to run this in gnu parallel ?

cat HaplotypeCaller.sh | parallel

#!/bin/bash

for i in *.bam; do
gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I $i -O $i.vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list; done

Answer: #!/bin/bash

do_one() {
gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I "$1" -O "$1".vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list
}
export -f do_one
# test that do_one works on a single bam file. Then:
parallel do_one ::: *.bam

If you do not want the output to be named `foo.bam.vcf.gz` but instead `foo.vcf.gz`:

#!/bin/bash

do_one() {
gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I "$1" -O "$2" --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list
}
export -f do_one
# test that do_one works on a single bam file:
# do_one foo.bam foo.vcf.gz
# Then:
parallel do_one {} {.} ::: *.bam
Question: Hi

I'm using plotHeatmap (library(scater)) to visualize the accuracy of cell type annotation output from SingleR. In total I got 28 labels, but as you can see on the top of the plot, the color different is hard to distinguish due to too many labels. I wonder whether there is a simple way to show only 5 to 7 labels in this plot, this is because I only care some specific cell types only, not all of them. For example "AT1 cells, AT2 cells and B cells" are the three cells I'm interested. How to only plot for only three labels?

Library(scater)
plotHeatmap(sce.data.combined.sct, order_columns_by = "labels", features=feature)

![heatmap output ][1]

[1]: /media/images/40053c91-45c3-473d-aa97-f50c8f7e

Answer: Gonna cheat and just recommend [dittoSeq](http://www.bioconductor.org/packages/release/bioc/html/dittoSeq.html). To do what you want will be something like:

```
dittoHeatmap(sce.data.combined.sct, genes = feature, order.by = "labels", annot.by = "labels",
cells.use = sce.data.combined.sct$labels %in% c("AT1 cells", "AT2 cells", "B cells"), scale = "none")
```

You may also be interested in setting `scaled.to.max = TRUE`, which is useful for single-cell data given their enrichment for 0s. Setting `scale = "row"` will plots z-scores instead. I also usually pass `complex = TRUE` and `use_raster = TRUE` with single cell data, as it makes the resulting image much easier to edit in Illustrator or such.
Question: Several pieces of software require this parameter.

Is counting the number of masked nucleotides in fasta files going to give a good approximate result?

If not, is there a simple way to do it?

Answer: EDIT: if you understand effective genome size as "mappable" genome size, than Devon is right, of course.

Assemblies only will provide you with a good size estimate if they are of really high quality. This is usually only the case for either model organisms or small, bacterial genomes.

Assemblies of larger genomes such as plant, animals etc., and in particular draft genomes usually do not contain a complete representation of a single haplogenome - which you would need to get your size estimate right. The reasons are that assembly algorithms usually cannot resolve all repeats and centromeric/telomeric regions, and also are prone to generate multiple sequences for different alleles of the the same region.

In my opinion, there are two better approaches:

1) Use the experimentally determined nuclear DNA (e.g. www.genomesize.com) content to calculate the haploid genome size. DNA content in pg can be directly converted into a bp estimate.

2) Use a k-mer based approaches to estimate the genome sizes form a high coverage NGS data set of your organism
Question: I am seeing strange behavior when comparing alignments generated by bwa mem 0.7.12 to two different reference genomes: hg19 and hs37d5 (basically hg19 plus additional decoy sequences). We have DNA-seq data. I noticed that when using hg19, there are very few alignments to the MHC region on chromosome 6, and even fewer of these have nonzero mapping quality. When using hs37d5, there are dramatically more alignments to the region and these have mostly high mapping quality scores. I have not observed this phenomenon anywhere else I've looked in the genome. The behavior is robust to multiple different choices of BWA parameters. Can anyone explain why the inclusion of the additional 35Mb of decoy sequences in hs37d5 would drastically improve the number and quality of alignments to this region of chr6?

Answer: I figured this out. It's because I was using the UCSC version of hg19 that includes some additional assemblies of MHC region haplotypes. Reads were preferentially mapping to those extra contigs instead of chr6.
Question: Hi everyone,

As my title states, I am looking for a database that has extensive collection of Arabidopsis structural RNAs. Your help will be much appreciated! THANKS!

Answer: [Ensembl Plants][1]. You can get *Arabidopsis t.**** ***gtf [here][2].

In this example you get CDS, exons, start/stop_codons

```
curl -s "ftp://ftp.ensemblgenomes.org/pub/release-24/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.24.gtf.gz" | gunzip -c |
cut -f 2 | sort | uniq -c

193 miRNA
1178 ncRNA
475382 protein_coding
1274 pseudogene
15 rRNA
71 snoRNA
13 snRNA
6277 transposable_element
697 tRNA
```

[1]: http://plants.ensembl.org/info/website/ftp/index.html
[2]: ftp://ftp.ensemblgenomes.org/pub/release-24/plants/gtf/arabidopsis_thaliana/