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: Theresa Scott's Intro deserves a mention although it can be a bit verbose.

Introduction to the
Fundamentals & Functionality of
the R Programming Language

For Bioconductor and short read analysis I would stick with Girke's cookbook which I have found amazingly helpful. I would gladly have paid $100 for access to it.

Sequence Analysis with R and Bioconductor

Have not read R in a Nutshell yet but from what I saw in Joseph Adler's Baseball Hacks and the initial reviews I have a good feeling about it.

I would avoid R Programming for Bioinformatics as that book is really geared toward package development.

I co-wrote an O'Reilly Short Cut, "Data Mashups in R", that is designed to be a little more fun than some of the scientific stuff out there while still exploring data manipulation in R, using packages, XML, web services, rudimentary plotting, and even some statistics. It costs $5 if you are not a Safari subscriber.

Data Mashups in R
Question: I can't find a clear definition that differentiates sequence assembly vs scaffolding.

From my understanding,

Assembly = joining reads into contigs

Scaffolding = joining contigs into scaffolds (using eg paired-end reads)

Does that sound right? It seems that assembly must be followed by scaffolding, but definitions of assembly don't even talk about scaffolding. Can you do assemble a whole genome with just "assembly"?

Answer: Contigs are sequences of overlapping (*contig*ous) reads. Paired-end (or mate-pair) reads can be used to determine the gap between two contigs. When you know the gap, you can make a scaffold, which is just the two contigs with Ns representing the gap in between.

EDIT: *Can you do assemble a whole genome with just "assembly"?*

Yes. Scaffolding won't give you more information about the actual bases anyways; it just tries to tell you how your contigs are ordered.
Question: Hello,

I have around 900 gene symbols and its Gene ID (HGNC) for which I need to extract all its corresponding OMIM phenotype. Is there a way to do it or a tool that does this ?

For example, *ABCA1* with HGNC ID 13666, the info to extract is listed here in column "[Phenotype][1]"

Thank you,

[1]: https://www.omim.org/entry/600046?search=abca1&highlight=abca1

Answer: I guess quering OMIM itself isn't possible. But you can use [BioMart][1] for it:

- Database: Ensembl Genes
- Dataset: Human Genes
- Copy&Paste your gene names to: Filters --> Gene --> Input external references ID list --> Gene Name(s)
- Select Attributes --> Gene --> Gene name and Phenotype --> Phenotype description and Phenotype --> Source name
- Click Results

You can export the result as CSV or XLS. You can then filter the result for "MIM Morbid" in the "Source Name" column to get these phenotype reported by OMIM.

fin swimmer

[1]: http://www.ensembl.org/biomart/martview
Question: Hi all,

In R, I have a data frame with physical position of some genes in bp on the left and intervals (again in bp) on the right. I want to check for each gene, whether it falls in the the interval of the same chromosome.

```
GeneID Chr Position Interval
x 1 18697251 1:3099740-8804450
y 1 19546617 1:3422930-8804450
z 1 3332236 2:2751757-4502486
a 2 3993537 2:6187995-8804450
b 2 9983384 2:4864334-18271005
c 2 11479025 2:8222941-18271005
```

For example: Is gene x on chr 1 comprised in one of the 2 intervals displayed for chr 1? I expect TRUE or FALSE.

I managed to do this giving the gene positions one by one, but the list is very big, I'd like this to be automated.

Any help much appreciated!

Answer: If you use `strsplit()`, you could test this directly. For instance:

> v <- '1:3099740-8804450'
> unlist(strsplit(v, '[:-]'))
[1] "1" "3099740" "8804450"

Then you can make a function out of it and print out the result of the desired tests via `apply()`:

> t <- read.table("some.table", header=T)
> interval <- function(str) unlist(strsplit(str, '[:-]'))
> apply(t, 1, function(x) (x[2] == interval(x[4])[1]) && (as.numeric(x[3]) > as.numeric(interval(x[4])[2])) && (as.numeric(x[3]) < as.numeric(interval(x[4])[3])))
[1] FALSE FALSE FALSE FALSE TRUE TRUE

We leave the chromosome name as a character vector to deal with cases like `X` and `Y`, and convert the position and interval position elements to numerics via `as.numeric()` to apply numerical relation operators.

The `apply` call can also be turned into a function:

> interval_test <- function(t) apply(t, 1, function(x) (x[2] == interval(x[4])[1]) && (as.numeric(x[3]) > as.numeric(interval(x[4])[2])) && (as.numeric(x[3]) < as.numeric(interval(x[4])[3])))

This function can be used to filter your table for rows that fall within the interval:

> t[interval_test(t),]
GeneID Chr Position Interval
5 b 2 9983384 2:4864334-18271005
6 c 2 11479025 2:8222941-18271005

Or, perhaps, used to filter for rows which do not fall within the interval, by using the `!` operator:

> t[!interval_test(t),]
GeneID Chr Position Interval
1 x 1 18697251 1:3099740-8804450
2 y 1 19546617 1:3422930-8804450
3 z 1 3332236 2:2751757-4502486
4 a 2 3993537 2:6187995-8804450
Question: I am trying to download and build the NCBI BLAST nr database, but I am running into the same problem constantly. The download seems to be successful, I don't receive any error messages and the size of the downloaded file is appropriate. Still, when I try to use the database in DIAMOND, I get the following error message:

> Error: Error reading input stream at line 360378522: Invalid character
> (>) in sequence

By running the following command, it seems like there is something wrong with the fasta.gz file:

zcat nr.gz | sed -n 360378500,360378540p

The problematic part looks like this:
>\>MBY6275551.1 hypothetical protein [Symbiobacterium thermophilum]
MRWSDVPVENKAQIAWGALVVVYLLGAQRPELIRPPVVAVFLLSAAAAFLELWLRSRGWPHLLWYDCIVWSALLTGMVVV
TGGRGSEVWAAYILMSLTAPVVLRRVAPYILLGVNVTAYGLIYLLYNPFGAPLDWGLLFLRIGTIFLVAYVVDRSTARER
QSHARAVALARSRVSELVQARDAERRRIAHDIHDWLGTGIIAPLRRLEIAARQSDVESCRRHVEEAADSLRRAHAELRRL
MENLHPHLLEQMGLAEALRAYLTDWGEEHGVAVHYHLTPGPEPPADAALALYRIQQEALNNCAKHADASQVWVTLELGAN
QVRLTLRDDGHGRPGRPGRGlacgTonARAVATIDPDPACSARQVGRWFYRAGKTGWQPLTRTSAGGCRDVHRVEPDFPEA
LCATAVRDILuFPkedLAELRRL
RQVGR>WJoIST TCP 71-13kei]NVFSIVGRWLPKQLGLKLAHCYEKNVEF guGIPH
MDCKCVa marETGNQGNSQGCEKV]OVseudal TNNRRPS
>GLDEF resYPKEGoPGTRLQVAAKATSYYPNQEMNILYSGGLDEERFPWTLNERVANDMRRmin1PNRAGRNRVANVVDAQLEETTPNISTAE

After this, I get the following error message:

>gzip: nr.gz: invalid compressed data--format violated

This is the same error message I get when I run

gunzip -t nr.gz

Interestingly I couldn't solve the problem by deleting and re-downloading the database and the position of the first invalid character varied between runs, so it doesn't seem like the problem is with the original file. I tested three ftp-downloading methods so far, I got the same end results with all three of them: FileZilla, wget and aria2. What am I doing wrong? What other methods can I try?
The commands I used to get the databases:

wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz

aria2c --continue=true --max-connection-per-server=4 --min-split-size=1M ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz

For FileZila, I simply navigated to the ftp site "ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz" and downloaded the nr.gz file

Answer: FOLLOW UP:

I finally solved this problem by downloading the file through FileZilla and setting "Transfer type" from "Automatic" to "Binary"
Question: Hi all!

I have hundreds of .csv files with the same dimensions and the same column 1 (gene_name). I would like to do an inner_join in the terminal (linux command line) in order to have a single matrix.

This is the .csv files `head SRR*_matrix_tpm.csv`
```
==> SRR5488793_matrix_tpm.csv <==
gene_name SRR5488793
NM_001042449.1 3752.38
NM_001043315.1 0
NM_001043316.2 1.31546
NM_001043358.1 0
NM_001043359.1 9.34679
NM_001043360.1 0.459982
NM_001043361.1 4.52898
NM_001043362.1 396.222
NM_001043363.1 101.929

==> SRR5488794_matrix_tpm.csv <==
gene_name SRR5488794
NM_001042449.1 5352.09
NM_001043315.1 0.0205977
NM_001043316.2 3.23229
NM_001043358.1 0
NM_001043359.1 15.0625
NM_001043360.1 0.833219
NM_001043361.1 11.0039
NM_001043362.1 147.638
NM_001043363.1 22.1844

==> SRR5570692_matrix_tpm.csv <==
gene_name SRR5570692
NM_001042449.1 5155.04
NM_001043315.1 1.74624
NM_001043316.2 18.6743
NM_001043358.1 0
NM_001043359.1 81.4916
NM_001043360.1 0
NM_001043361.1 16.1719
NM_001043362.1 34391.8
NM_001043363.1 3315.22
```
I would like to have something like this:

```
==> matrix.csv <==
gene_name	SRR5488793	SRR5488794	SRR5570692
NM_001042449.1	3752.38	5352.09	5155.04
NM_001043315.1	0,0	0.0205977	174624,0
NM_001043316.2	131546,0	323229,0	186743,0
NM_001043358.1	0,0	0,0	0,0
NM_001043359.1	934679,0	150625,0	814916,0
NM_001043360.1	0.459982	0.833219	0,0
NM_001043361.1	452898,0	110039,0	161719,0
NM_001043362.1	396222,0	147638,0	34391.8
NM_001043363.1	101929,0	221844,0	3315.22
```
`Many thanks!`

Answer: Unix `join` command can only be used to join 2 tables together at a time. You will need a few piped commands to get the joined result.

```
$ join -1 1 -2 1 SRR5488793_matrix_tpm.csv SRR5488794_matrix_tpm.csv | tr ' ' '\t' | join -1 1 -2 1 - SRR5570692_matrix_tpm.csv | tr ' ' '\t' > joined.tsv`
$ cat joined.tsv
gene_name	SRR5488793	SRR5488794	SRR5570692
NM_001042449.1	3752.38	5352.09	5155.04
NM_001043315.1	0	0.0205977	1.74624
NM_001043316.2	1.31546	3.23229	18.6743
NM_001043358.1	0	0	0
NM_001043359.1	9.34679	15.0625	81.4916
NM_001043360.1	0.459982	0.833219	0
NM_001043361.1	4.52898	11.0039	16.1719
NM_001043362.1	396.222	147.638	34391.8
NM_001043363.1	101.929	22.1844	3315.22
```
*EDIT: the above only works if the files are sorted for the 'key', in this case column 1*

If you have hundreds of files with the same length and order, you can use the `cut` command to create

```
cut -f1 SRR5488793_matrix_tpm.csv > rownames.txt
for CSV in *csv ; do
cut -f2 $CSV > $CSV.cut
done
paste rownames.txt *cut > joined.tsv
```

Alternatively in R there is a `join_all()` function in the `plyr` package.
Question: I am using freebayes and bwa-mem together to align and call variants on my sequencing reads. But I am doing some editing to the reads before aligning and calling variants. I want some loci to be ignored based on my edits.

So if for example I have a base whose WT is G, but for a given read I see an A, if I change that A to an N will the sequence then be properly aligned to the genome but freebayes will not report any info for that locus, such that I won't see the variant, but I also won't see a read supporting the WT sequence?

Answer: If your sequencing data has quality scores (FASTQ), you should set the quality of the base to a very low value (depending on how your quality scores are encoded). The mapper will know to ignore it.

You can set the base to N, or you can even leave it as-is and hope the mapper maps it with mis-matches (which it should), but dropping the quality is probably the best way.

Even better would be if you edited the genome you're mapping too ;)
Question: Hi

I'm using STAR to perform alignment and my goal is to do the DEG analysis in the future. The parameter I set up as follows:

Step 1:
```
STAR --runThreadN 12 \
--runMode genomeGenerate \
--genomeDir genomedir \
--genomeFastaFiles ./ref/GRCm39.primary_assembly.genome.fa \
--sjdbGTFfile ./ref/gencode.vM27.primary_assembly.annotation.gtf \
--sjdbOverhang 100
```
Step 2:
```
STAR --runThreadN 12 \
--runMode alignReads \
--twopassMode None \
--genomeDir genomedir \
--sjdbGTFfile ref/gencode.vM27.primary_assembly.annotation.gtf \
--sjdbOverhang 100 \
--readFilesIn "R1_001_trim.fastq.gz" "R2_001_trim.fastq.gz" \
--readFilesCommand zcat \
--outFileNamePrefix "./align_output/try1_" \
--outSAMtype BAM SortedByCoordinate \
--quantMode GeneCounts
```
But the output of genecount file is a bit weird:
```
N_unmapped	1657081	1657081	1657081
N_multimapping	8053375	8053375	8053375
N_noFeature	12418438	43724858	13355695
N_ambiguous	1629488	39643	430809
ENSMUSG00000102693.2	0	0	0
ENSMUSG00000064842.3	0	0	0
ENSMUSG00000051951.6	49	2	47
ENSMUSG00000102851.2	0	0	0
ENSMUSG00000103377.2	10	0	10
ENSMUSG00000104017.2	13	0	13
ENSMUSG00000103025.2	4	1	3
ENSMUSG00000089699.2	0	0	0
ENSMUSG00000103201.2	2	0	2
ENSMUSG00000103147.2	1	1	0
ENSMUSG00000103161.2	7	0	7
ENSMUSG00000102331.2	20	0	20
ENSMUSG00000102348.2	0	0	0
```

Why the gene index at the first column has ".2" or ".3" in the end? like "ENSMUSG00000102693.2", Shouldn't be just "ENSMUSG00000102693"? I use reference from Gencode by the way

Answer: That's the version number of that ensembl gene ID. When a gene model is changed to an appreciable extent the version number increases.
Question: I'am trying to open index.html from JBrowse after running ./setup.sh.
I have just loaded sample_data of volvox.
But some of the tracks are not loading file.
Showing "Error : HTTP 200 when fetching file"

Please help me in this
Thank You.!

Screenshot : https://photos.app.goo.gl/NQsa5SJC5w7APfHk8

Answer: Thank you everyone..

Actually i wanted to implement browse in Flask framework, and by just adding "JBrowse" folder in static directory of my app solved this problem.

Thanks again for your valuable responses.
Question: Hi,

Dear Community,

I am not sure to post this type of question here, but I observed some posts related to setting up Bioinformatics lab.

I need to ask one thing. Do you know any company or charity NGO or organization that donates old / refurbished computer servers (having good configuration for Bioinformatics work) to universities? As people donate old scientific equipment, is this happened in Bioinformatics field?

Regards,
Waqas

Answer: You can try various local NGO's or international companies involved in electronic donation projects.

Though, it will be hard to get high end functional computers which could support heavy bioinformatics tasks (image processing, 3D modelling, clusters etc.) but these organizations provide functional computers which are used for under-privileged people/children in developing and under-developed countries for basic learning and communication purposes.

Few sources, I'd like to mention :

[Computer Aid][1] - Versatile selection

[World Computer Exchange][2]

[Interconnection][3] - Subsidized computer store for low income/under-developed.

Another way is to get sponsored from big companies like Novartis, Roche etc. Most of them they have a sperate division where they run off-shoot projects caring for environment, health care and other sectors. One can write them a proposal and sometimes, they can sponsor the whole project.

Good Luck!

[1]: http://www.computeraid.org/apply.html
[2]: http://www.worldcomputerexchange.org/apply-computers
[3]: http://interconnection.org/lowincome/desktops.html
Question: I have a protein motif or site, which I like to identify in an DNA sequence (multiple fasta file). The motif is N-X-S/T (X!=P), which means Asn, followed by any amino acid but not Pro, followed by Ser or Thr. Also X should not be STOP. So I would like to find all the 3 codon combinations for this site in DNA (9 nucleotides).

I was first thinking of getting the motif written in DNA using IUPAC coding, but that seemed not possible. Writing out all possibilities seems like a too hard task, so I thought there might be a tool which can do this? Any suggestions?

Answer: Haven't tried it, but you could do a 2-tiered grep-approach. Make sure your fasta is not line-wrapped.

grep -o "AA[CT][ACGT]\{3\}\([AU]C[ACGT]|AG[CU]\)" fasta_file | grep -v "[ACGT]\{3\}CC[ACGT][ACGT]\{3\}".

Assuming `Asn = AAY = AA[CU]`, `Ser = UCN, AGY = UC[ACGT], AG[CU]`, `Thr = ACN = AC[ACGT]`, and `Pro = CCN = CC[ACGT]`,
the first part should match all peptides `N-X(traditional = all amino acids)-S/T`, the second should get rid off the ones that contain proline in the central position. I am not sure about whether you have to use an additional set of `\(\)` in the first expression.
Question: I am trying to perform a GSEA for ten biochemical properties to check which properties are enriched for my set of genes. I am using the Windows application (v4.2.2) with GSEAPreranked package. However, when I perform the analysis some of them returns a p-value of "0.0". When I checked the user guide it states:

"In the GSEA report, a p value of zero (0.0) indicates an actual p value of less than 1/number-of-permutations. For example, if the analysis performed 100 permutations, a reported p value of 0.0 indicates an actual p value of less than 0.01. For a more accurate p value, increase the number of permutations performed by the analysis. Typically, you will want to perform 1000 permutations (phenotype or gene_set). (If you attempt to perform significantly more than 1000 permutations, GSEA may run out of memory.)"

I tried increasing the permutation to 10,000 and 100,000 which still yielded a "0.0" value. Is there a way to calculate the absolute p-value which I can use to later calculate FDR q-value?

Answer: You may want to consider switching to [fgsea](http://bioconductor.org/packages/release/bioc/html/fgsea.html) which has algorithmic improvements to GSEA for better p-value estimation [[PREPRINT](https://www.biorxiv.org/content/10.1101/060012v3)]. It helps to resolve most cases where the p-value is really small.
Question: This is my samtools flagstat output for the alignment of a paired-end sample using HISAT (filtered with MQ<1):

```
45084184 + 0 in total (QC-passed reads + QC-failed reads)
4717987 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
45084184 + 0 mapped (100.00%:-nan%)
40366197 + 0 paired in sequencing
20273012 + 0 read1
20093185 + 0 read2
39146254 + 0 properly paired (96.98%:-nan%)
39644363 + 0 with itself and mate mapped
721834 + 0 singletons (1.79%:-nan%)
260722 + 0 with mate mapped to a different chr
235361 + 0 with mate mapped to a different chr (mapQ>=5)
(END)
```

There's quite a few questions about flagstat on Biostars, but I haven't seen the secondary, supplementary or duplicates flag before. Assuming that this is HISAT specific, are is secondary the number of mapped reads that get mapped using the second (non-global) HISAT index? What would supplementary mean? And are duplicates the amount of reads that map non-uniquely?

Answer: These are defined more formally in the SAM specification, though perhaps the wording there isn't great.

Secondary: One of the many places a multimapper can align. Note that multimappers will have one primary and 0 or more such secondary alignments.

Supplementary: For chimeric/fusion/non-linear alignments, this is the location of one part of the alignment. For reads that align in a chimeric fashion, one segment will be designated as primary and the remainder supplementary.

Duplicate: If you've marked possible PCR duplicates, then this will be set. The definition of a duplicate is somewhat dependent on the tool used (N.B., aligners don't typically set this flag, it's down by picard's markDuplicates command or similar).
Question: I am trying to pull regions from a bunch of BAM files on an online server. I'd like to pull the reads mapping to a certain 1kb or so chunk and download them for analysis. They are far too massive to download them all, and it's impractical even to wget them one at a time and pull the regions out using samtools (tried it, and it worked, but it took forever). Since I'll have to do this for a number of regions that I won't know in advance, I need a better way.

I noticed that samtools is capable of running 'samtools view' off of a web address. Sadly, this data is protected behind an https server, which samtools doesn't know how to handle. I notice that IGV is able to read the BAM files of the net by asking for my login and querying specific regions only that I bring up to view, but I don't have a way of automating the process on hundreds of files.

Does anyone have any ideas of how to run something like samtools view on specific regions over an https connection?

Answer: samtools, and other direct-BAM-access programs like IGV, are capable of opening local files as well as remotely and publicly served files. these remote locations include ftp and http protocols, but unfortunately do not include encrypted transfer protocols such as scp, ssh nor https. it's not a matter of authentication, which is solved on http and ftp, but of handling data encryption which is far more complicated. the only way you may work with all that BAM files you are interested in is either asking the server managers to open them through http, or either downloading them all and dealing with them locally.
Question: Dear all,

I need to perform a large BLAST search and I am using blastn in remote from the terminal. However, this takes way too long to give an answer and I have been thinking of creating a local database to speed the analysis. How can I download the all nr/nt repository? I see there is one [here][1] for the RefSeq. Would be this good? would it be already indexed or shall I create the index with makeblastdb?

Thank you

[1]: https://www.ncbi.nlm.nih.gov/refseq/

Answer: You can use

wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz

wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz

These are fasta files, they are not indexed. You should use the `makeblastdb` command to index that.

You might also want to browse ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA and check what other databases are available.
Question: Hi guys, I have a question about how to define a gene to be expressed in RNA Seq analysis.
Is it better to use log2(rpkm) or rpkm without log2 transformation to define a gene to be expressed? I know that the log is only a transformation so that log2(rpkm) will let to consider as expressed a gene with a higher rpkm with respect to the simple rpkm. However I think that rpkm alone leads to consider as expressed genes that have too few reads.

Some ideas about?

Thank you in advance

Answer: If you just want to filter out genes of low expression prior to performing downstream analyses, then just go by the RPKM values and eliminate genes with mean RPKM < 10 or < 20, and / or those genes with a high frequency of 0 or NA values.

If your aim, however, is to define, for example, an expression signature for a particular tissue or condition, then transform your RPKM dataset to Z-scores via *zFPKM* package (R Programming Language) and then set cut-offs:

- Z > 2 = expressed
- Z < -2 = not expressed

Those cut-offs can, of course, be modified.

Note that neither RPKM, FPKM, nor CPM data is suitable for differential expression analysis, and neither are the logged transformations of these.

Kevin