cannin/biostars_qa huggingface dataset:
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: SimpleR
Question: Dear Biostars, Hi

I want to run **blastx** against **NCBI nr** for annotating **thousands** of *de novo* transcripts. I want to do it as fast as possible.

I have heard that one way is to use Amazon cloud.
So, would you kindly help and introduce me any step-by-step manual or resource in this regard?

(e.g: where to begin? quantity of resources I need to rent from cloud and related costs? should I download the Nr database locally or downloading it into Amazon cloud? or uploading my transcriptome assembly into the Amazon? is there any pre-installed and ready to use Cloud for Bioinfirmatics (some BioCloud)? . . . and so on).

Thanks

**Notes:**

0- I have checked some [Biostars posts][1] but could not find any simple guidance.

1- I have no computational resources at home (no server and ...)

2- I do not want to use *Blast2GO Pro* CloudBlast (why? it is expensive for me!)

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

Answer: Have you looked at the [NCBI BLAST on Cloud help][1]? There is step-by-step guidance available using links at the left.

That said, if you were looking to do this for thousands of transcripts then you should use [DIAMOND instead][2]. Regular blast+ may take too long and cost more.

There is also [PAUDA][3].

[1]: https://ncbi.github.io/blast-cloud/
[2]: https://github.com/bbuchfink/diamond
[3]: https://ab.inf.uni-tuebingen.de/software/pauda
Question: Hello!

I am currently analyzing a set of RNA-seq data. I first used StringTie to calculate the TPM values of reads that are mapped to the genome (a non-model organism). Then, for the ~5% of reads that were not mapped to the genome, I used Trinity to *de novo* assemble the reads and then used Kallisto to calculate the TPM of these transcripts.

Are the TPM values calculated by StringTie and Kallisto comparable? For example, can I combine the two lists of TPM values and look at the top 30 expressed transcripts?

Thank you very much!

Answer: No, I don't think concatenating your TPMs like this is a good idea: both data sets are normalized to 1million, however the number of reads is very different between them.

If you want to do something like this, I would:

1) map reads to genome

2) used unmapped reads to assemble transcriptome

3) extract / get transcripts from genome

4) concatenate assembled transcriptome and extracted transcripts into new reference transcriptome

5) use kallisto to estimate TPMs / counts for new reference transcriptome
Question: Hi all,

I am using htseq-count for getting the read count of my genes from RNASeq data.While I have got the read counts but I am getting the warning below:

```
Warning: 35062120 reads with missing mate encountered.
35138401 SAM alignment pairs processed.
```

I am guessing this could be because of how I sorted my sam file.I did not specify any parameters while sorting.This was the command I used(gnu parallel). {1} is used to specify my input files.

samtools sort - -m 40G {1/.}.sorted

My question is do I have to sort by any specific parameter such as -n ?For example do I need to rerun sorting and then produce results of htseq-count.

Thanks in advance

Ron

Answer: First check how many reads have proper mates. May be using samtools flagstat

If you find many reads with missing mates, like shown in htseq warning, there are no issues with htseq. If you find most of the reads have their mates aligned, then

Sort the samfile either by position ( default ) or name ( -n ). Then use latest version of htseq and set the -r option to eithe name ( -r name [default] ) or position ( -r pos ) based on how you have sorted the bam file. Make sure you have used the strandedness (-s) option appropriately.
Question: Hi, apologies for this basic questions, I new in NGS quality control. I have been check my NGS data (Illumina - HiSeq 2500 2*100pb) using FastQC after trimming Nextera Adaptater with bbduck (BBTool package) using trimming overlap (ktrim=r k=25 mink=11 hdist=1 tpe tbo). And the checks fails Per base sequence content, Per sequence GC content and Kmer Content. My question is should I be trying remove the first base (~15) of my sequence? When I try this, Kmer Content fails again but at the end of the sequence. Should I remove again at the end?
After my goal is to assemble them with a De Bruijn assembler and thus i'm afraid of over-represent kmer.

Thanks for the help.

![enter image description here][1]
![enter image description here][2]

[1]: https://i.imgsafe.org/498873a.png
[2]: https://i.imgsafe.org/8d597c0.png

Answer: This [post][1] from Dr. Simon Andrews (author of FastQC) should be the last word about this question that keeps coming up in one form or other.

[1]: https://sequencing.qcfail.com/articles/positional-sequence-bias-in-random-primed-libraries/
Question: Hi everyone,
I would like to get SNP positions of a particular strain of mouse (JF1) from VCF file deposited in ftp server of Sanger Institute. The link of the file is as follows: ftp://ftp-mouse.sanger.ac.uk/REL-1807-SNPs_Indels/mgp.v6.merged.norm.snp.indels.sfiltered.vcf.gz.
It would be very helpful if somebody can suggest me a tool or way to do this.
Thanks
Ankit

Answer: wget -q -O - "ftp://ftp-mouse.sanger.ac.uk/REL-1807-SNPs_Indels/mgp.v6.merged.norm.snp.indels.sfiltered.vcf.gz" |\
bcftools view --samples "JF1_MsJ" |\
bcftools query -i 'INFO/AC>0' -f '%CHROM %POS\n'
Question: I am in the process of updating some annotations using the Ensembl Biomart.

Using the biomaRt package in R I do the following:

mart = useMart("ensembl", dataset="rnorvegicus_gene_ensembl")
transcript_mart = getBM(
attributes = c(
"ensembl_transcript_id",
"ensembl_gene_id"),
filters="ensembl_gene_id",
values="ENSRNOG00000047860",
mart=mart)

The result is empty. I check on the website and get a result:

Ensembl Gene ID Ensembl Transcript ID
ENSRNOG00000047860 ENSRNOT00000073930

So the thing is: The website uses Ensembl Genes 71 (the example gene was added with the new rat genome assembly in 70).
I check the biomaRt package and it is the current 2.16.0.

listMarts()
biomart version
1 ensembl ENSEMBL GENES 69 (SANGER UK)

tells me though that biomaRt 69 is current.

listDatasets(mart)
dataset description version
12 rnorvegicus_gene_ensembl Rattus norvegicus genes (RGSC3.4) RGSC3.4

And uses the old genome assembly.

If I check the [Biomart registry](http://www.biomart.org/biomart/martservice?type=registry)

Still the old version 69 marts are listed.

So the question is: Am I missing something or is the biomaRt package/registry simply out of date?

Answer: Yes, everything that you state is correct and it seems that version 69 is still the Bioconductor/biomaRt default.

You can use release 71 as explained in this post from the BioC mailing list: [Useful information about Ensembl release 71 mart databases](http://grokbase.com/t/r/bioconductor/134bx70fe3/bioc-useful-information-about-ensembl-release-71-mart-databases).

mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
host = "www.ensembl.org",
path = "/biomart/martservice",
dataset = "rnorvegicus_gene_ensembl")

transcript_mart <- getBM(attributes = c("ensembl_transcript_id", "ensembl_gene_id"),
filters = "ensembl_gene_id",
values = "ENSRNOG00000047860", mart = mart)

transcript_mart
# ensembl_transcript_id ensembl_gene_id
# 1 ENSRNOT00000073930 ENSRNOG00000047860
Question: The Bam format was developed at a time when I/O (disk / network) was not a limiting factor for genomics data analysis. Currently it is. Especially when you have pipelines like GATK that process BAM files in multiple steps and read and write bam files for every step (instead of in memory streaming).

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

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

Answer: Neither of the solutions you have is a full alternative: HDF and CRAM are generic binary formats that could store any data. Neither is a specification (nor is it obvious that a SAM file transformed into HDF would actually be more performant than a BAM format). The reduced BAM format is not a generic format either, it is a simplified format that targets one specific use case.
Question: Hello,

I have a dataset with a variety of samples that vary in Age (Young/Old) and Sex (M/F).

I'm interested in testing a few hypotheses, including (Q1) "What genes are DE as a product of Sex?" and (Q2) "What genes are DE as a result of the interaction of Age and Sex?"

To answer Q1, I originally imported data like so:

# Attempt 1
myData <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable_forAllMySamples,
directory = pathToHTSeq,
design = ~Sex)
dds <- DESeq(myData)

This produced a very large DE gene list.

Later, I redid this analysis with a different design matrix including interaction and contrasts, like so:

# Attempt 2
myData <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable_forAllMySamples,
directory = pathToHTSeq,
design = ~Sex+Age+Sex:Age)
dds <- DESeq(myData )
res <- results(dds, contrasts=c('Sex', 'M', 'F'))

However, this produced a MUCH smaller list of DE genes.

My understanding is that pulling out the contrasts should look for the main effect of Sex in my dataset (so, Sex effects regardless of Timepoint).

I had expected that would be the same as if I just made design matrix ~Sex, but it looks like that isn't the case. **Why is that?**

It that because Attempt 2's design matrix "controls" for Age and any interaction effects, but Attempt 1 does not? Can anyone help me understand a bit better what is being tested in Attempt 1, or point me towards resources to strengthen my understanding of what that was doing?

Possibly relevant: When I PCA plotted my rlog-normalized data, the data clustered very well by Sex, and less well by Age.

Thank you very much for your help!

Answer: Things like slightly imbalanced group sizes (in this case, the numbers of males and females at each age) as well as the difference power increasing with sample size are the prime causes for this. I should note that the results with just `~sex` as the design likely have more false-positives, since they're not accounting for the confounder of age (as you astutely surmised).

It's pretty common for samples to cluster strongly by sex, its effect isn't as variable as something like age.
Question: Hi,

I am learning STAR using two sets of 1000 RNA-seq reads and mapping to the human genome.

But am having a problem: When the STAR mapping run finishes no Alignemnt.out.sam file is created.

I am using the following:

/home/programs/STAR-STAR_2.4.2a/bin/Linux_x86_64_static/STAR \
--runMode alignReads \
--limitGenomeGenerateRAM 31000000000 \
--limitIObufferSize 150000000 \
--genomeDir /home/test_index \
--runThreadN 8 \
--outFilterMismatchNmax 2 \
--outFilterMultimapNmax 10 \
--genomeLoad NoSharedMemory \
--alignIntronMin 21 \
--outStd SAM \
--outSAMattributes Standard \
--outSJfilterReads Unique \
--sjdbGTFfile /home/gene_models/Homo_sapiens.GRCh38.76.gtf \
--readFilesIn /home/reads/short_reads_human.fastq /home/ddowling/lsf_test/short_reads_human_2.fastq

STAR finishes and creates the following files in the working directory:

_STARgenome
SJ.out.tab
Log.std.out
Log.progress.out
Log.out
Log.final.out

The 'Alignment.out.sam' is missing.

The Log.std.out reads:

14:21:23 ..... Started STAR run
14:21:23 ..... Loading genome
14:22:21 ..... Processing annotations GTF
14:22:49 ..... Inserting junctions into the genome indices
14:25:03 ..... Started mapping
14:25:13 ..... Finished successfully

So I am unsure where the problem is.

Thank you.

Answer: You used --outStd SAM, which would have outputted the Aligned.out.bam file to your stdout.

From [the docs][1]:

> --outStd default: Log string: which output will be directed to stdout (standard out) SAM: alignments in SAM format (which normally are
> output to Aligned.out.sam file), normal standard output will go into
> Log.std.out

So check the other log files for a SAM file, and if it's not there, read the manual next time! (and be glad you didn't make the same mistake with a much slower aligner ;)

EDIT: But upvote for your question because you asked it excellently :)
[1]: https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf
Question: [note: cross-posted on GATK forum - still awaiting a response]

I have a GVCF (generated using GATK's HaplotypeCaller w/ -ERC GVCF parameter) of 36 related samples and would like to determine the (potentially *de novo*) variants that are unique to each sample. Short of creating 36 N-1 GVCFs for discordance testing, or extracting individual sample VCFs for subtraction, is there a straightforward method to obtain the desired information?

Answer: merge your g.vcf files with CombineGVCFs https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_CombineGVCFs.php

and then I would use bioalcidaejdk to find the sample containing a singleton: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

java -jar dist/bioalcidaejdk.jar -e 'stream().forEach(V->{final List L=V.getGenotypes().stream().filter(G->G.isHet() || G.isHomVar()).collect(Collectors.toList());if(L.size()!=1) return;final Genotype g=L.get(0);println(V.getContig()+" "+V.getStart()+" "+V.getReference()+" "+g.getSampleName()+" "+g.getAlleles());});' input.vcf

* `stream():` get a stream of variants
* `final List L=V.getGenotypes().stream() :` get a stream of genotpes for variant V
* `filter(G->G.isHet() || G.isHomVar()).` get only genotypes HET or HOM_VAR
* `collect(Collectors.toList());` convert to list
* `if(L.size()!=1) return;` list must contain a singleton
* `final Genotype g=L.get(0);` the singleton
* `println(V.getContig()+" "+V.getStart()+" "+V.getReference()+" "+g.getSampleName()+" "+g.getAlleles()); ` print the singleton
Question: Dear all,

Now we used ANNOVAR to annotate our SNP detected. But the result of ANNOVAR is the single-letter Amino acid abbreviation. It's like this:

The results I want is like this:

So does anyone knows that which annotation software can feed back this 3-letter amino acid abbreviation? Thanks in advance.

Best wishes,

Answer: update:

Actually snpEff can feed back 3-letter Amino acid abbreviation. And the results are quite good.
Question: Hello,

I'm trying to merge peak files in a way that I can't find so far.

I have 16 peak files concatenated in my_peaks.bed

I've tried to use bedtools merge following this example:

$ cat my_peaks.bed
chr1 100 300
chr1 200 500
chr1 0 700

$ bedtools merge -i my_peaks.bed -c 1 -o count
chr1 0 700 3

But actually the kind of output that I want would be the following:

$ bedtools merge -i my_peaks.bed -c 1 -o count
chr1 0 100 1
chr1 100 200 2
chr1 200 300 3
chr1 300 500 2
chr1 500 700 1

I've looked into bedops option --partition which would give :

$ bedops --partition my_peaks.bed
chr1 0 100
chr1 100 200
chr1 200 300
chr1 300 500
chr1 500 700

But then I don't know how to keep the peak count.

Thanks for help

Answer: You're very close; just pipe the partioned elements to BEDOPS *[bedmap][1]* to count peaks over partitions:

$ sort-bed my_peaks.unsorted.bed > my_peaks.bed
$ bedops --partition my_peaks.bed | bedmap --echo --count --delim '\t' - my_peaks.bed
chr1	0	100	1
chr1	100	200	2
chr1	200	300	3
chr1	300	500	2
chr1	500	700	1

Note that you first need to sort your peaks. Using (BEDOPS) `sort-bed` will be faster than Unix `sort` and will sort on the stop column, where there are matches on the start column, which `sortBed` does not do. You only need to sort once, at the start.

[1]: http://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html
Question: In the Cuffdiff outputs I can see the FPKM of each gene in the groups, but how can I get the FPKM of each gene in each one of the samples in the groups???

Answer: As juafonso_bio has already found out, the FPKM values in the output of Cuffdiff can be found in the file genes.read_group_tracking.
The replicates are identified by numbers, starting at 0.
The correspondance between the replicate number and the sample name can be found in the file read_groups.info.

Here is an example.

-bash-4.1$ more genes.read_group_tracking
tracking_id	condition	replicate	raw_frags	internal_scaled_frags	external_scaled_frags	FPKM	effective_length	status
ENSMUSG00000000001	shCas1	2	8025	6474.79	6474.79	56.6616	-	OK
ENSMUSG00000000001	shCas1	0	6508	6724.17	6724.17	58.844	-	OK
ENSMUSG00000000001	shCas1	1	4740	5932.42	5932.42	51.9153	-	OK
ENSMUSG00000000001	shRnf2_388	0	7579	7162.23	7162.23	62.6775	-	OK
ENSMUSG00000000001	shRnf2_388	1	6783	7142.23	7142.23	62.5025	-	OK
ENSMUSG00000000001	shRnf2_388	2	7982	6855.53	6855.53	59.9935	-	OK
ENSMUSG00000000001	Casz1v2	0	7719	8440.91	8440.91	73.8673	-	OK
ENSMUSG00000000001	Casz1v2	1	6487	7135.16	7135.16	62.4406	-	OK
ENSMUSG00000000001	Casz1v2	2	7712	7013.6	7013.6	61.3768	-	OK
ENSMUSG00000000001	pSiren	1	6116	6064.53	6064.53	53.0714	-	OK
ENSMUSG00000000001	pSiren	0	7557	6877.29	6877.29	60.1839	-	OK
ENSMUSG00000000001	pSiren	2	6151	6762.49	6762.49	59.1793	-	OK

-bash-4.1$ more read_groups.info
file	condition	replicate_num	total_mass	norm_mass	internal_scale	external_scale
../../results/tophat/shCas1_1/accepted_hits.bam	shCas1	0	3.24875e+07	3.5031e+07	0.967852	1
../../results/tophat/shCas1_2/accepted_hits.bam	shCas1	1	2.96395e+07	3.5031e+07	0.798999	1
../../results/tophat/shCas1_3/accepted_hits.bam	shCas1	2	4.18947e+07	3.5031e+07	1.23942	1
../../results/tophat/shRnf2_388_1/accepted_hits.bam	shRnf2_388	0	3.64417e+07	3.5031e+07	1.05819	1
../../results/tophat/shRnf2_388_2/accepted_hits.bam	shRnf2_388	1	3.28269e+07	3.5031e+07	0.949703	1
../../results/tophat/shRnf2_388_3/accepted_hits.bam	shRnf2_388	2	4.05007e+07	3.5031e+07	1.16432	1
../../results/tophat/Casz1v2_1/accepted_hits.bam	Casz1v2	0	3.22445e+07	3.5031e+07	0.914475	1
../../results/tophat/Casz1v2_2/accepted_hits.bam	Casz1v2	1	3.26844e+07	3.5031e+07	0.909159	1
../../results/tophat/Casz1v2_3/accepted_hits.bam	Casz1v2	2	3.91855e+07	3.5031e+07	1.09958	1
../../results/tophat/pSiren_1/accepted_hits.bam	pSiren	0	3.79887e+07	3.5031e+07	1.09883	1
../../results/tophat/pSiren_2/accepted_hits.bam	pSiren	1	3.63825e+07	3.5031e+07	1.00849	1
../../results/tophat/pSiren_3/accepted_hits.bam	pSiren	2	3.15436e+07	3.5031e+07	0.909576	1
Question: I am using BWA to align RNA-seq data from a miRNA experiment.

I want to compare myself to another article - so I am using the same techniques its authors used.

That is using fastq-mcf to clip the adapters, and BWA (Burrows Wheeler Aligner) to align the reads.

The reads are not aligned to the entire genome. Instead they are aligned to the mature miRNA forms, downloaded as fasta files from miRBase (example below).

After doing the alignment I get a SAM file (example below).

I want to get the count of each miRNA from that SAM file.

The way I thought of doing the count, is to count the number of rows that map to each miRNA in the SAM file. I would only take the rows in which FLAG (second column) & 4 == 0. That is - the read is mapped. The miRNA identity will be computed from the RNAME (third column).

Is this the right way to do the counting?

To begin with, I thought of using htseq-count, but the SAM file doesn't really contain the coordinates of the miRNAs in the reference sequence, so htseq cannot work.

Samples lines from the fasta file of the mature miRNA forms:

```
>mmu-let-7g-5p MIMAT0000121 Mus musculus let-7g-5p
TGAGGTAGTAGTTTGTACAGTT
ACTGTACAGGCCACTGCCTTGC
>mmu-let-7g-3p MIMAT0004519 Mus musculus let-7g-3p
```

Sample lines from the SAM assembly file:

```
HISEQ:47:C6FUWANXX:1:1101:2080:1992 4 * 0 0 * * 0 0 TTGGCAATGGTAGAACTCACACCGAAA

Answer: You probably want to ignore unmapped and secondary alignments:

samtools view -F 260 foo.bam | cut -f 3 | sort | uniq -c | awk '{printf("%s\t%s\n", $2, $1)}' > counts.txt

or something along those lines.
Question: if genomic DNA is extracted from tissues, and the goal is to find variants such as SNP, indels, CNVs-which is the preferred method- whole exome sequencing or whole genome sequencing and why?

Under what cases can we use RNA sequencing to find variants?

Answer: > if genomic DNA is extracted from tissues, and the goal is to find variants such as SNP, indels, CNVs-which is the preferred method- whole exome sequencing or whole genome sequencing and why?

SNPs and short indels can be called both using WES and WGS. However, you can only call variants if they're sequenced, which with WES is limited to the exome, obviously. If you are interested in non-coding variants you need WGS. CNVs and other SVs are hard with short-read WGS and nearly impossible with WES.

But of course there is also a difference in price, so your budget will play a role as well.

> Under what cases can we use RNA sequencing to find variants?

You should never design an RNA-seq experiment with the aim to find variants. However, if RNA-seq is what you have you can look at variants, but take into account that it's far from optimal and won't give you the best results.
Question: Hello.

I am looking for a download from the EggNOG FTP site for the most current mapping of orthologous groups (NOGs) to their Gene Ontology (GO) terms.

[I see where I can query for the Go terms if I know the NOG number][1], but I am working in the other direction: I have GO terms and want to know to which NOGs they correspond.

Additionally, I see where [mapper][2] is useful if I want to get GO terms for a previously unmapped group of sequences. However, I am wanting to work with NOGs in the EggNOG, so those should already be mapped and mapping them again seems unnecessary and inefficient.

Help is appreciated! And thanks.

[1]: http://eggnogdb.embl.de/#/app/api
[2]: http://eggnogdb.embl.de/#/app/emapper

Answer: After several hours of downloading and checking through files from the EggNOG FTP, I finally found this one:

http://eggnogdb.embl.de/download/eggnog_4.5/all_OG_annotations.tsv.gz

The file just under 1G and seems to contain mappings to other annotation type databases.

One line of data contains as follows and is available in a compressed tab-delimited file.

> 0RWQ4 homNOG 4 RNA binding motif protein 43 [u'B']
> {"Molecular Function":[["GO:0097159","organic cyclic compound
> binding","IEA",4,"100","100"],["GO:0000166","nucleotide
> binding","IEA",4,"100","100"],["GO:0036094","small molecule
> binding","IEA",4,"100","100"],["GO:0005488","binding","IEA",4,"100","100"],["GO:1901363","heterocyclic
> compound binding","IEA",4,"100","100"],["GO:1901265","nucleoside
> phosphate binding","IEA",4,"100","100"],["GO:0003676","nucleic acid
> binding","IEA",3,"75","98.4"],["GO:0003723","RNA
> binding","IEA",1,"25","92.5"]]} []
> {"SMART":[["RRM",3,"75","98.4"]]}
> 9601.ENSPPYP00000014337,9606.ENSP00000331211,9593.ENSGGOP00000008646,9598.ENSPTRP00000021421

The "0RWQ4" is the EggNOG unique identifier for the protein OG. As we might predict, not all OGs have GO mappings.