cannin/biostars_qa huggingface dataset:
Question: Hi all,

I am looking for a piece of advice related to a efficient and accurate way of obtaining a multiple sequence alignment of all the 3'UTRs of a genome. I'm aiming to detect conserved and repetitive patterns that could be targets for regulatory elements, but I am encountering memory limitations (classic multiple alignment tools -such as Muscle- cannot handle files that contain more than 500 sequences or their size is over 10 Mb).
I would hugely appreciate any ideas or feedback.
Thank you very much in advance!

Answer: Personally I don't think that trying to align all 3'UTR of genes in a genome will result in anything decent! There is way too much variation (== too little conservation) to result in a sensible multiple alignment.

Moreover, simply doing a global alignment is probably technically not even feasible. Lengths of the 3'UTR are too variable and if you're looking for motifs those are too small (insignificant) to be visible in a global multiple alignment.

A better approach (= what most people will do ) is to first group genes (UTRs) in biological significant clusters (eg. co-expression, same pathway, ...) and then try to analyze the UTRs in those smaller groups. Even so, the simple MSA approach will then not even give much I assume. You're better trying software specifically for that purpose (= motif detection), I'm thinking of phylogenetic footprinting, phylogenetic shadowing, motif sampling, RSA-tools ....
Question: I've noticed poor overlap between number of differentially expressed genes between CLC's 'Empirical Analysis of DGE' (edgeR test) versus Cuffdiff. CLC maps 20% more reads on average, and finds 200-400% more differentially expressed genes. It is to my understanding that CLC uses a proprietary read mapper that will usually give a higher percentage, while Tophat's probabilistic model spends more time resolving overlapping ends to reduce random placement. However, it is less obvious how to relax stringency in Tophat's parameters to account for the high amount of variation in the hybrids we sequenced versus our reference genome. Any insight would be greatly appreciated.

Answer: You'll probably want to experiment with different parameters, but Tophat is extremely flexible with regards to its parameters. Here are some you'll want to check out. From the [manual][1] (which is well worth your time to read through to get an idea of the knobs you can manipulate):

`-N` maximum number of mismatches a read can have without being discarded

`-a` How many bases must be mapped to each side of a junction

`-m` How many mismatches are allowed in an anchor region of a spliced alignment

`--report-secondary-alignments` Don't report only the best alignments

`-g` The maximum total alignments that will be reported for a read

`--segment-length` You can reduce this from the default of 25 to cut the reads up into smaller chunks, which would make the mapping more generous

`--segment-mismatches` How many mismatches a given segment is allowed to have

You can also try using `--b2-very-sensitive` for a quick setup

Without knowing more about your experiment I can't recommend specific starting parameters, but it's worth varying the parameters of interest by large amounts and comparing the results side-by-side in [IGV][2] or [IGB][3].

[1]: http://ccb.jhu.edu/software/tophat/manual.shtml
[2]: http://www.broadinstitute.org/igv/
[3]: http://bioviz.org/igb/index.html
Question: Hi, I am just curious. Can we use STAR aligner to align the paired-end reads as a single-end? What are the consequences of the output of such alignment?

Answer: Yes, STAR can align paired-end reads as single-end, just run STAR two times, one with `--readFilesIn
R1.fastq`, another with `--readFilesIn R2.fastq`.

> What are the consequences of the output of such alignment?

First and foremost, the mapped reads will be considered as single-end. In addition, the proportion of multi-mappers may increase.
Question: Hi All,

Supposed I have a .gtf file with the exons of all genes in a given genome. I would like to extract the first exon coordinate of each gene from the .gtf file. How should I get started?

Regards,
Ziyi

Answer: As your question is about getting started, I would suggest to look for code snippets or libraries that parse GTF files to get an idea about handling a GTF file.

For example:

http://www-huber.embl.de/users/anders/HTSeq/doc/tour.html#tour

https://github.com/ctokheim/PrimerSeq/blob/master/gtf.py

But a quick and dirty way would be:

curl https://raw.githubusercontent.com/roryk/DEXSeq/master/inst/python_scripts/dexseq_prepare_annotation.py | python - genes.gtf out.tmp​
grep "exonic_part_number \"001\"" out.tmp | less -S

This gives all the first exonic parts of a gene, assuming a standard gtf file format.
Question: We are often suggested to sort the input bed file by "sort -k1,1 -k2,2n" in order to invokes a memory-efficient algorithm designed for large files, for example, bedtools intersect ( http://bedtools.readthedocs.org/en/latest/content/tools/intersect.html)

But this is slow for a large file with >20G in size. Any quick-sorting problem around?

Here is a solution I can think of:

inputbed=$1
awk -v inputbed=$inputbed '{print $0 >> inputbed"."$1}' $inputbed
for i in $inputbed.*; do sort -k2,2n $i > $i.s & done
# when it's done
sort -m -k1,1 -k2,2n $inputbed.*s > $inputbed.sorted

The basic idea is to split the large file into small files by the first key, then sort each of them by the second key, and finally use "sort -m" to merge the sorted ones. It can save time by parallel sorting individual small file. But you need a way to track when the individual sorting is done.

I can also read the records into a large hash table, and then sort the key before output. Below is an example of sorting by 1st column (using a hash of hash can sort by two keys):

perl -e 'while (<>) {$l=$_; @a=split("\t", $l); push(@{$HoA{$a[0]}}, $l);}{foreach $i (sort keys %HoA) {print join("", @{$HoA{$i}});}}'

But again, this method requires large memory.

I am wondering if any quick-sorting program around that you guys can recommend. Thanks

Edit1: Use the parallel sorting from GNU coreutils with "sort --parallel=N" (change the number of sorts run concurrently to N) option. Also set the main buffer size and traditional locale.

LC_ALL=C sort --parallel=24 --buffer-size=5G -k1,1 -k2,2n input.bed > sorted.bed

Edit2: Use the sort-bed tool from BEDOPS

sort-bed --max-mem 5G input.bed > sorted.bed

Answer: You might not need to go through the pain of rolling your own.

There are multithreaded sort commands.

have a look @ stackoverflow

http://stackoverflow.com/questions/6365176/multi-threaded-sort
Question: Hi there,
I am trying to analyze 10X genomics output with Seurat. After reading in data using "Read10X" function and create an Seurat object, I am stuck at this step:

> mito.genes <- grep(pattern = "^MT-", x = rownames(x = my_data@data), value = TRUE)
For the data I have, row names of the expression matrix are gene names ("Sox17" "Mrpl15" "Lypla1" "Tcea1" "Rgs20" etc) or ensembl ID, where no genes begin with "MT-".

Does anyone know the mitochondria gene list Seurat uses for quality control? Or I should generate a mitochondria gene list on my own? I am working on mouse samples.

Many thanks for your help!

Answer: Looks like you have a mouse genome. Instead of `^MT-`, you can look for `^Mt-` or just add `ignore.case = TRUE` to the `grep`.
Question: Hello,

I wonder what is the best way/package to extract ~3k bps regions surrounding a certain given location in R. E.g. the location can be Human Chromosome 4 Position 50000. Then I'd like to have the sequence from chr4:47000-53000 of hg19/hg38.

Thank you.

Answer: You can use [biomaRt][1]

GENES = useMart("ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl", host = "useast.ensembl.org")

YouRange<- getBM(attributes = c("chromosome_name", "start_position", "end_position"), filters = c("chromosome_name","start","end"), values = list(chromosome_name = "4", start = 50000-3000, end = 50000 + 3000), mart = GENES)

[1]: https://bioconductor.org/packages/release/bioc/html/biomaRt.html

You can use `getSequence` get sequences.

getSequence(chromosome = 4, start = 50000-3000, end = 50000+3000,type='hgnc_symbol',seqType = 'gene_exon_intron', mart = GENES)

`getBM` helps you get the annotation information from Biomart, `getSequence` helps you get the sequence you need.

This is a simple example for example question, you can check details of biomaRt to solve your specific question.
Question: How can I input reads from a file descriptor rather than an actual text file? I have created four small test FASTQ files to identify the problem; two R1 files and two R2 files. If I use process substitution, Clumpify fails. If I firstly combine R1 reads into a text file on disk and R2 reads into another file and use them, then Clumpify works. The reason I am exploring this is because I have a data set of 51 whole genome DNA samples. 47 of them were sequenced in one lane and four of them were sequenced in two lanes. The command I used is

clumpify.sh -Xmx1g t=1 in1=<(cat ${R1[@]}) in2=<(cat ${R2[@]}) out1=$R1output out2=$R2output dedupe subs=1

This results in the error

Starting cris 0.
Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: 1

The only difference for the command which works successfully is `in1=merged1.fastq` and `in2=merged2.fastq`

Done!
Time: 0.971 seconds.
Reads Processed: 4000 4.12k reads/sec
Bases Processed: 600k 0.62m bases/sec

I got this idea from the STAR RNA-seq aligner, which permits `--readFilesIn <(gunzip -c reads.fastq.gz)`, for example.

On the university's computing cluster, version 37.98 of bbmap is installed. I could use the `/scratch/` directory if all else fails.

Answer: I think the problem is `clumpify.sh` won't work with process substitution, regardless of the arrays. Here is a solution similar to https://www.biostars.org/u/25721/ but using arrays, and interleaving with `paste` and `awk` taken shamelessly from this [gist][1]:

paste <(zcat ${r1[@]}) <(zcat ${r2[@]}) \
| paste - - - - \
| awk -v OFS="\n" -v FS="\t" '{print($1,$3,$5,$7,$2,$4,$6,$8)}' \
| clumpify.sh -Xmx1g t=1 int=t in=stdin.fastq out1=$R1output out2=$R2output dedupe subs=1

[1]: https://gist.github.com/nathanhaigh/4544979
Question: Hi all,

I have been using R to create and manipulate phylogenetic trees for a while, but I was wondering what the differences are between the rooting functions in different packages, `ape::root`, `ggtree::reroot` and `phytools::reroot`.

I have attached an example. Every tree looks slightly different, and the tree created by `ggtree::reroot` seems best to me. Could anyone explain the differences?

tree <- read.tree(text = "(t6:0.3806255851,(t7:0.5440872659,(t5:0.6203179485,t1:0.7089423968):0.3421749519):0.2618428892,((t4:0.5886056426,t3:0.8614832051):0.09094934678,(t8:0.6766210056,t2:0.1412485428):0.4451906278):0.2356684168);")

par(mfcol=c(2,2), mar=c(0.5,0.5,1.3,0.5))

plot(tree)
nodelabels()
title("original")

tree1 <- ape::root(tree, node = 11)
plot(tree1)
nodelabels()
title("ape::root")

tree2 <- ggtree::reroot(tree, 11)
plot(tree2)
nodelabels()
title("ggtree::reroot")

tree3 <- phytools::reroot(tree, 11)
plot(tree3)
nodelabels()
title("phytools::reroot")

![Example][1]

Any explanation would be greatly appreciated.

Best wishes,

Michael

[1]: https://i.imgur.com/WM5XPgf.png

Answer: Just to clarify, the reason the tree from ggtree 'looks' better is because it is still bifurcating - but a bifurcating tree rooted at an internal node should contain a trifurcation at the root (as do the trees from phytools & ape). To obtain a bifurcating tree from rerooting an existing bifurcating tree we must position the root along an edge - but an edge is defined by two nodes, not one. Furthermore, we also need to specify where on the edge the rerooting must take place. In phytools this can be done using the argument *position*, or by selecting it interactively by setting interactive=TRUE (in which case we can click on the position along the edge where we would like the tree to be rooted).
Question: Is there a Picard Tool to filter one or more reads through the following :

FLAG (second column)

Mapping Quality (MAPQ 5th Column)

TLEN (9th Column)

??

Answer: Hypothetical scenario: Suppose you'd want to get all reads with flag 83=paired reads|first in pair|read reverse mapped| proper pair reads with mapping quality > 30 from a bam file and return output to a bam file:

samtools view -h input.bam | awk 'substr($0,1,1) == "@" || $2 == 83 || $5 >= 30 {print}' | samtools view -bS - > output.bam

How it works?:

samtools view -h input.bam - decodes bam and prints sam file with header to screen
awk 'substr($0,1,1) == "@" - keep all header lines
|| $2 == 83 || $5 >= 30 {print} - the actual condition you would want to check. for unmapped reads you could use GAWK with bitwise operation and($2, 4) == 4 || and($2, 16) == 16, or something like that.
samtools view -bS - pipe input S=sam to output b=bam and the "-" mediates samtools to take input from the screen
> output.bam - simple unix redirection to output file

Basically, you can use the awk line to filter anything you want.
Question: Hello all,

I wanted to get a poll of what distros are currently popular with researchers working in bioinformatics and related biological data mining fields.

Mainly:

What linux distribution do you currently (or have previously) use(d)?

What software and language (packs) do you use on a daily basis? -

Where does this distro outperform the rest? Where does it fall short?

I'll start; I currently use Ubuntu and have for the past two years. I mainly use vim with BioPerl scripts. Ubuntu, which is debian based makes installing foreign libraries almost painless. It falls short when it comes to software updates.

Answer: Currently in use: Arch Linux

Software packages:

R (statistical analysis)

Scala (programming tasks)

Perl / BioPerl (simple pipelines)

blast / bwa (local alignments)

LaTex / Inkscape (presentation and visualization)

Processing (Visualization with java)

Pros and cons: (decide yourself what is pro and what is con)

Very clean and lightweight

Hard to configure

Very configurable

Best documented linux

Your learn how linux works

Great community

Always current versions of packages

Sliding releases

Conclusion:

Perfect distro if you know what you are doing.
Question: In PEER-paper they write:

> When RNA-seq estimates are used for transcript abundance, we recommend using DESeq to estimate library sizes and variance-stabilize expression data sets.

And I got confused by DESeq2 and DESeq. In DESeq (varianceStabili... depends on the estimateDispersions):

```r
cdsBlind = estimateDispersions( cds, method="blind" )
vsd = varianceStabilizingTransformation( cdsBlind )
```

In DESeq2:

vsd <- varianceStabilizingTransformation(dds)

I am getting the same results from `varianceStabilizingTransformation()` if I use DESeq functions as `DESeq()`, `estimateSizeFactors()`, `estimateDispersions()` before or not.

```r
dds = DESeqDataSetFromMatrix(countData=histone_m, colData=Design, design=~condition)
pik = DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

vst=varianceStabilizingTransformation(pik)
head(assay(vst))
WEN1 WEN2 WNN1 WNN2
1 8.568387 8.206238 7.843942 8.086705
2 8.962147 8.688945 8.320699 8.750248
3 8.135942 7.859804 7.843942 7.878520
4 7.321410 6.863988 7.165757 6.908507
5 6.474284 6.379227 6.420262 6.351544
6 6.561254 6.468225 6.555806 6.861555

vst=varianceStabilizingTransformation(dds)
head(assay(vst))
WEN1 WEN2 WNN1 WNN2
1 8.568387 8.206238 7.843942 8.086705
2 8.962147 8.688945 8.320699 8.750248
3 8.135942 7.859804 7.843942 7.878520
4 7.321410 6.863988 7.165757 6.908507
5 6.474284 6.379227 6.420262 6.351544
6 6.561254 6.468225 6.555806 6.861555

cds=estimateSizeFactors(dds)
cds = estimateDispersions( cds)
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates

vst=varianceStabilizingTransformation(cds)
head(assay(vst))
WEN1 WEN2 WNN1 WNN2
1 8.568387 8.206238 7.843942 8.086705
2 8.962147 8.688945 8.320699 8.750248
3 8.135942 7.859804 7.843942 7.878520
4 7.321410 6.863988 7.165757 6.908507
5 6.474284 6.379227 6.420262 6.351544
6 6.561254 6.468225 6.555806 6.861555
```

So, it seems it does not matter if I use DESeq function before or not. Does it make sense? Also, once I estimate SizeFactors and varianceStabilizingTransformation, how should I plug it into PEER?

Answer: From my understanding, there is negligible difference between the Variance Stabilizing Transformation functions of DESeq and DESeq2. I would suggest you use DESeq2 because it is the most up-to-date version.

```
library(peer)

# assign vst data to variable called expr_set
expr_set = assay(vst)

# the PEER tutorial recommends k == 25% of "number of samples"
k = ncol(expr_set) * .25

# create empty peer model
model = PEER()

# set number of "hidden factors" searched for
PEER_setNk(model,k)

# PEER ask NxG matrix, where N=samples and G=genes
PEER_setPhenoMean(model,as.matrix(t(expr_set)))

PEER_setAdd_mean(model, TRUE)

# These can be changed to the tolerance you desire
PEER_setTolerance(model, .001)
PEER_setVarTolerance(model, 0.0005)

# usually doesn't go over 200 iterations
PEER_setNmax_iterations(model, 1000)

# run peer
PEER_update(model)
```
Question: Of course, as a bioinformatician, I am aware of many large-scale open-source bioinformatics datasets, such as

- The ENCODE consortium (www.encodeproject.org, RNA-Seq, ChIP-Seq and so on),
- The Roadmap Epigenomics consortium (www.roadmapepigenomics.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The IHEC consortium (www.ihec-epigenomes.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The TCGA/ICGC consortia (www.cancergenome.nih.gov, www.icgc.org, large-scale cancer data, DNA-Seq, RNA-Seq, etc.) and
- The LINCS consortium (www.lincscloud.org/l1000 , gene expression for more than a million of different perurbation experiments).

I am wondering, however, what other wonderful datasets, the are both large and open-source, are currently available. That might include things like RNA-Seq, Chip-Seq, Bisulfite-Seq, whole genome sequencing, WGAS, and many other assays (not necessarily NGS-related, though that is what I am mostly looking for).

Also things like the (neural) connectome of certain species (in any event large data) could be of interest.

There are quite some GEO datasets that at least partially fulfill these requirements, but most are simply having to few data samples in order to be interesting to me.

Your comments are greatly appreciated!

Answer: SRA maybe? http://www.ncbi.nlm.nih.gov/sra

2.2 Petabases open source, 1.8 Petabases authorized access. (btw, TCGA is authorized access not open access)
Question: There is a lot of tools out there, very useful for command line usage, and very widely used in Bioinformatics, which rapidly turns out to be annoying (may be, sometimes) if we are writing a pipeline that cares about I/O connection and each tasks' exit status.

I am writing a pipeline using samtools, and samtools turns out to be a little bit annoying in the I/O management, because sometimes it generates an output file, but you don't really explicitly name that file. Sometimes, other tools don't even prompt for output files, or some other tools ask users to provide paths literally which adds up more turnarounds that need to be introduced and this can be a bit frustrating.
Here is an example using samtools. I am wrapping a call to samtools on a file that does not exist, the command is failing but the exit code is still zero which is a bit misleading if we care about reporting the status of the entire pipeline, which means here it means that the sort went Ok and this will trigger other tasks, which is wrong

I am cross posting to CodersCrowd as well with a code you can run on the browser: http://coderscrowd.com/app/codes/view/288. You can see that either the status coming back from the docker image and the one coming back from the python interpreter itself (which is basically samtools exit status) is being zero, and it shouldn't be.

![< image not found >][1]

[1]: http://s12.postimg.org/izkv3yeu5/Screen_Shot_2014_12_02_at_2_18_11_PM.png

Answer: Thank you all for your contributions.

To all those who will be interested in this kind of problems, the overall message that should be retained from this discussion is : when there is a problem in the source, it should be treated in the source, whatever wrapper, codes, tricks that you may end up doing, will most likely be not reproducible at all.

I end up cross posting to samtools mailing list as well and here is the final thought :

```
You are testing an old version of samtools. In the particular case of samtools sort's exit status
$ samtools-1.0 sort S222.bam S22.sorted; echo "exit code $?"
[E::hts_open] fail to open file 'S222.bam'
[bam_sort_core] fail to open file S222.bam
exit code 1
```

So as you can see, it was a minor bug, and that was fixed in the last version of samtools

Thank you all for contributing to this discussion
Question: Hello, I'm searching for a way to calculate the `ESTIMATE tumor purity,` using genes expression matrix.
I found something close, https://github.com/mwsill/RFpurify

However this doesn't seem suitable for RNA-seq data, is was trained on methylation array data.

Any suggestions for such a package in R?

Answer: Just in case anyone faces this problem in the future, here is the link for the updated `estimate package` on github:

https://github.com/KaiAragaki/tidyestimate

for some reason this package is hidden. After days of searching I have finally found it.
Question: Is it possible to pass paired end data on fastqc commandline?

The

./fastqc --help

doesn't mention how to do it. It is possible to pass multiple files like

./fastq file1 file2 file3

but it doesn't seem to differentiate somehow SE with PE.

On [another post][1] I read that FASTQC is not suited for paired end. However I got confused because inside the **Configuration/contaminants.txt** file the reads are paired end.

[1]: https://www.biostars.org/p/65001/

Answer: Normally fastqc will analyse the files for the forward and reverse reads separately.