cannin/biostars_qa huggingface dataset:
Question: I need to convert `blast` output in tabular format (`format 6`) to `bed` format. I'm interested only in `qstart` and `qend` coordinates (the 7th & 8th columns respectively). As far i know `bed` format is 0-based, but `blast` probably uses 1-based coordinates. Am i right?

So to do so with `AWK` i need to execute the code below:

cat balst_tabular_output.txt | awk '{print($2"\t"$7-1"\t"$8)}' > blast_otput.bed

Answer: Yes, BLAST is 1-based and your conversion should be right (keeping qend 1-based makes it compatible with the open interval style that bed uses.)
Question: Hello,

I wanted to obtain only the "mapped" reads as an output of the STAR.
I forgot to delete "**--outSAMunmapped Within** " and all of my mapped output files also contain all unmapped reads,too.
Data size is huge and re-mapping properly wiil take so much time...

**How can I fix the BAM files which were supposed to be mapped reads but also include unmapped reads in it ?**

#!/bin/bash
mkdir /mnt/data/Toxo_scan/GBR_Male/ToxoMap
while read -r line
do
mkdir ToxoMap/$line
echo $line" -> Running STAR - Toxo now"
STAR --runThreadN 12 --alignIntronMax 1 --outSAMunmapped Within --outSAMtype BAM SortedByCoordinate --genomeDir "/mnt/data/Toxo_scan/toxo_genome" --readFilesIn "/mnt/data/Toxo_scan/GBR_Male/sickle/"$line"_1_clean.fastq" "/mnt/data/Toxo_scan/GBR_Male/sickle/"$line"_2_clean.fastq" --outFileNamePrefix "/mnt/data/Toxo_scan/GBR_Male/ToxoMap/"$line"_" --outReadsUnmapped Fastx
done

Answer: To get your mapped reads:

samtools view -b -F 4 alignment.bam > mapped.bam

To get your unmapped reads:

samtools view -b -f 4 alignment.bam > unmapped.bam

The `-f` filters based on bitwise samflags. Flag 4 (or `0x0004`) is for unmapped reads. `-F` does the inverse.
Question: hi,

i have a correlation matrix of genes but i am going to remove alphabet from rows and columns to prepare my file for matlb i did like below

> mycounts <- read.table("tmatnorm_rld500.txt", header = T, sep = "\t", row.names=1)
> head(mycounts[,1:3])
AT1G01060 AT1G01170 AT1G01180
AT1G01060 1.0000000 0.3193360 0.6272994
AT1G01170 0.3193360 1.0000000 0.3178880
AT1G01180 0.6272994 0.3178880 1.0000000
AT1G01260 0.2658624 0.3588030 0.2542234
AT1G01380 0.6178751 0.1561297 0.7052692
AT1G01490 0.5990499 0.6129301 0.6424225
> rownames(mycounts) <- NULL
> head(mycounts[,1:3])
AT1G01060 AT1G01170 AT1G01180
[1,] 1.0000000 0.3193360 0.6272994
[2,] 0.3193360 1.0000000 0.3178880
[3,] 0.6272994 0.3178880 1.0000000
[4,] 0.2658624 0.3588030 0.2542234
[5,] 0.6178751 0.1561297 0.7052692
[6,] 0.5990499 0.6129301 0.6424225
> colnames(mycounts) <- NULL
>
> head(mycounts[,1:3])
[,1] [,2] [,3]
[1,] 1.0000000 0.3193360 0.6272994
[2,] 0.3193360 1.0000000 0.3178880
[3,] 0.6272994 0.3178880 1.0000000
[4,] 0.2658624 0.3588030 0.2542234
[5,] 0.6178751 0.1561297 0.7052692
[6,] 0.5990499 0.6129301 0.6424225

> write.table(mycounts, file = "correlation.txt", dec = ".", sep = "\t", quote = FALSE, row.names=F)
> mycounts <- read.table("correlation.txt", header = T, sep = "\t")
> head(mycounts[,1:3])
V1 V2 V3 <---this line
1 1.0000000 0.3193360 0.6272994
2 0.3193360 1.0000000 0.3178880
3 0.6272994 0.3178880 1.0000000
4 0.2658624 0.3588030 0.2542234
5 0.6178751 0.1561297 0.7052692
6 0.5990499 0.6129301 0.6424225

but as you consider, i have something in colnames?? then how I can get rid of alphabet and only have the similarity values please?

thank you

Answer: Add `col.names=FALSE` to the `write.table()` command and you should have a matrix with no row or column names.
Question: Hi guys,

I have a multi-fasta like this

>Citrobacter_freundii_D8_6645..17576
gtgatcgtcaagaaggttaagaacccgcagaaggcagca
>Enterobacter_hormaechei_35012_3830..23574
atggacgatagagaaagaggcttagcatttttatttgcaatt

And I would like to eliminate the numbers flanking .., to have an output like this

>Citrobacter_freundii_D8
gtgatcgtcaagaaggttaagaacccgcagaaggcagca
>Enterobacter_hormaechei_35012
atggacgatagagaaagaggcttagcatttttatttgcaatt

Since the number are variable, I guess just creating a command to remove x characters from the end of the fasta-header won't be enough.
Thanks!

Answer: If the example is representative, then you basically intend to keep the first three elements that are separated by `_`. If so, do:

awk ' $1 ~ /^>/ { split($0,a,"_"); print a[1]"_"a[2]"_"a[3];next} {print}'

Command splits every line that starts with `>` at the `_` and then simply prints the first three separated by `_` again. Obviously that only works if all fasta headers look like the ones you showed.
Question: Using R I want to generate a list of genes only (without the accompanying text) from a pathway.

For example. If the input pathway is [KEGG prostate cancer](https://www.genome.jp/dbget-bin/www_bget?hsa05215),
I want my output to be a .csv list of the genes in that pathway. I tried:

library("KEGGREST")

keggGet("hsa05215")[[1]]$GENE

but that gives me a list of the number and gene description along with the gene symbol and I want a list consisting of the gene symbol alone.

How do I get this?

Thank you.

Answer: I ended up solving it myself. Hopefully this helps anyone who has a similar need.

library("KEGGREST")

#Get the list of numbers, gene symbols and gene description
names <- keggGet("hsa05215")[[1]]$GENE
#Delete the gene number by deleting every other line
namesodd <- names[seq(0,length(names),2)]
#Create a substring deleting everything after the ; on each line (this deletes the gene description).
namestrue <- gsub("\\;.*","",namesodd)
#export the vector as a csv
write.csv(namestrue, file = "hsa05215",quote = F, row.names = F)
Question: Hello all

I want remove some snp from map and also genotype from ped file. I used this command in plink:

plink --file FILE --exclude SNP-LIST-FILE

I think the some SNP in map file should be deleted and also genotype in ped, but when I check the number of raw in map and the column of ped file its same as before exclude

is there any one know my snp in snplist excluded or not?

thanks

Answer: Your command is correct, but you should add two options: `--out` and `--recode`. If you don't specify these, plink executes the commands and prints a summary to the screen (and to the .log file), but does not actually write the results to a ped file nor overwrites your original data.

The `--recode` option specifies that you want to write a .ped file as output. The `--out` option specifies the suffix of the output file:

plink --file mydata --exclude snp_list --recode --out filter_data

Also, check that the the SNP-LIST-FILE is formatted correctly. It must contain the ids of the SNPs that you want to remove.
Question: Does O, P, or Q stands for different things ?

I have looked at this webpage, but it didn't answer the quesion
https://www.uniprot.org/help/accession_numbers

Answer: Confirmation from UniProt: Indeed accession numbers have no meaning and just serve as unique identifiers.

This is the short answer, and below is some background:

Originally, accession numbers started with P, later with O, P and Q,
but due to the large increase in the number of protein sequences in
UniProtKB, we had to extend the existing accession number format by
allowing the first character to be any of the 26 letters (instead of
only O, P and Q). To avoid assigning accession numbers identical to
those which have been used by the International Nucleotide Sequence
Database, the extension in the first position went along with a
restriction in the third position which can only be a letter. See
https://www.uniprot.org/news/2006/11/14/release

The most recent extension of the UniProtKB accession number format
took place in 2014 to
allow for 10-character accession numbers:
http://www.uniprot.org/help/2014/06/11/release

The "new-style" accession numbers do not have any different meaning at
all, except that they implicitly convey the information that the entry
they identify has been created since June 2014 - however the reverse
is not necessarily true, i.e. not all new entries have 10-character
accessions. In any case, it is best not to interpret ACs in any way,
but just consider them as unique identifiers for a UniProtKB record.
Question: Hi every one i am doing some plamid genome assembly with spades. After assembly i used SSPACE for scaffolding. But there are some gap each of the draft genome. I can fill the gap by PCR. But i want to reduce the number of gap insilico?? Can anyone suggest how to reduce the gap ?? If i do mapping the reads with contigs will it be give any promising result???

Answer: There are a few in silico tools that can fill in the gaps within scaffolds;

GapFiller (developed by myself, I'm also the developer of SSPACE):
http://www.baseclear.com/landingpages/basetools-a-wide-range-of-bioinformatics-solutions/gapfiller/
http://genomebiology.com/2012/13/6/R56

IMAGE2:
http://sourceforge.net/apps/mediawiki/image2/index.php?title=Main_Page

GapCloser (part of SOAP):
http://sourceforge.net/projects/soapdenovo2/files/GapCloser/

Regards,
Boetsie
Question: I tested two different options while running HTSeq-Count, `-s no` and `-s reverse`. This are the results:

For `-s no`:
> __no_feature	435592

> __ambiguous	953159

> __too_low_aQual	0

> __not_aligned	0

> __alignment_not_unique	8164048

For `-s reverse`:
> __no_feature	573728

> __ambiguous	410510

> __too_low_aQual	0

> __not_aligned	0

> __alignment_not_unique	8164048

For the option `-s reverse` there are lower ambiguous values but higher no_feature than for `-s no`.
As far as I know this option depends on the construction of the library, but when they gave me this sequences they didn't mention it. All I know is that it's was constructed under a Illumina protocol and that it was a Paired-End experiment of RNA-Seq from peach (Prunus persica).
I'm inclined to think that less ambiguous values are just better, even than with more no_feature values.

So, which one it's right?

===============

Edit:
This are the results from `-s yes`:
> __no_feature	41467373

> __ambiguous	506

> __too_low_aQual	0

> __not_aligned	0

> __alignment_not_unique	8164048

Answer: This is a pretty common method to determine the strandedness of a library and is essentially what RSeQC is doing. In this case, `-s reverse` is the correct setting (it's also the majority of what's produced these days). The general reasoning is:

1. If it's an unstranded library, then each of the stranded methods will have ~2x more `_no_feature` counts than the `-s no` setting.
2. If it's a stranded library, one of the stranded setting will have "slightly" higher `_no_feature` counts (because reality is annoying like that) and the other will have vastly higher `_no_feature` counts.
Question: I am confused about an approach to evaluate enrichment and was wondering if you could help me understand if what I am doing makes sense.

I have counts from regions that overlap histone markers in a dataset and I would like to know if histone markers are enriched in my dataset compared to a random set using a randomization approach.

I have created 1000 similar datasets and created counts for number of regions overlapping histones in these null datasets. In some cases, this distribution is normal, but in others it is not.

In cases where the distribution of the proportions from the null datasets is **normal**:

I can use these data to find the mean, standard deviation, and degrees of freedom and compare this distribution (mean, sd and df) to my observed count. Is this correct?

```r
# "sim.null" is a normal distribution of counts of overlaps that I get from 1000 simulations (this is matched on my original dataset for some features)

sim.null= rnorm(sd=0.001, mean=0.01, n=1000)

# I would like to compare it with the counts I get from my dataset

observed = 0.0125

t = mean(sim.null)-observed / (sd(sim.null)/sqrt(1000))

2*pt(-abs(t),df=999)

# Is this the same as doing this?

t.test(sim.null, mu=observed, alternative="two.sided")$p.value
```

In cases where the null dataset is **not normal**: would maybe a Fisher exact test be appropriate?

Thank you very much, any suggestions are very appreciated!

Answer: What you meant to write was:

Tstat = (mean(sim.null)-0.0125)/(sd(sim.null)/sqrt(1000))

The p-value will then essentially match (the t.test will tell you <2e-16, which is the smallest value it'll output).

For non-normal distributions, either use the appropriate distribution, or just use an empirical distribution (the latter is required for messy distributions):

sum(sim.null>=0.0125)/1000

Using theoretical distributions (e.g., the normal distribution) will give much nicer p-values if they're appropriate for your dataset.
Question: I am looking at the ENCODE ATAC-seq pipeline: https://www.encodeproject.org/pipelines/ENCPL035XIO/

They have two different steps:

- "call nuclease accessible regions using FSeq" (in PDF) or "open chromatin region identification" (on diagram)
- "call nuclease accessible peaks using Homer" (in PDF) or "peak calling" (on diagram)

Regardless of the tool used, what is the difference between "regions" and "peaks"? I would think those are the same thing (in this context, a set of loci where the reads accumulate).

Answer: I received a very helpful clarification after emailing ENCODE directly:

> Nuclease accessible regions tend to be long, e.g. 10 kb or longer.
> This was clear even in the early papers on DNase sensitivity
> (mid-to-late 1970's; Groudine and Weintraub). These accessible regions
> can contain entire genes or even clusters of genes. Within the
> nuclease accessible regions, some localized DNA segments are so
> readily cleaved that double-strand breaks are generated at that
> position in a substantial fraction of the cells in the population.
> These are the DNase-hypersensitive sites (DHSs) first mapped by Carl
> Wu (late 1970's). I see the Fseq "regions" as the equivalent of
> nuclease accessible regions, and the Homer "peals" as the equivalent
> of DHSs.
>
> If you look at the signal track for DNase-seq or ATAC-seq, you see
> broad regions of signal that are significantly above the background.
> Within those regions, you see localized peaks, often many peaks per
> region. Fseq calls the broad regions, and we use Homer to call the
> localized peaks. MACs can be used for peak calling as well, Anshul
> Kundaje is doing that. You can see similar analyses in the work from
> John Stamatoyannopoulos for DNase-seq. I think Hotspots are like
> regions, and DHSs are peaks confined to a defined length.
Question: Hi guys I am trying to learn use R for RNA-SEQ analysis and I am trying some functions and guides from bioconductor, however I am getting some errors.

http://www.sthda.com/english/wiki/rna-seq-differential-expression-work-flow-using-deseq2#introduction

One of them I quickly found a answer, because they change the name of the fuunction 8(.

But this error I didn't find any answer still:

txdb <- makeTranscriptDbFromBiomart( biomart="ENSEMBL_MART_ENSEMBL",
dataset="hsapiens_gene_ensembl" )
Error in makeTranscriptDbFromBiomart(biomart = "ENSEMBL_MART_ENSEMBL", :
could not find function "makeTranscriptDbFromBiomart"

I import the library as due

library( "GenomicFeatures" )
But got that error.
Any tip?
thanks

Answer: I got it done with makeTxDbFromBiomart...
well I guess they change the functions names a lot in R packages.

Thsnks
Question: I want to query the molecular weight of multiple proteins from R given uniprot ids. I thought of using biomaRt or UniProt.ws, but didn't see which field I should query.

There are online tools performing the calculation given a uniprot ID (PIR, ExPASy), but you have to manually enter the input and save the output.

I want the MW that appears in UniProt when looking in "Sequences" under Mass (for example for the entry http://www.uniprot.org/uniprot/Q9Z2Y8)

Answer: On the uniprot site you can ask for the predicted mass column.

Click on the columns button. click on the check-box labeled mass under the sequence category.

Save, search for your accession using the syntax below and get tab delimited output back.

http://www.uniprot.org/uniprot/?format=tab&query=accession:Q6GZX4&columns=id,entry+name,mass
Question: Q: What should I put in the read groups in my BAM files?

Yes, I've read about the "read groups" thing on biostars (e.g. Picard provides a very useful tool to add/replace read groups). But I think I've missed something very fundamental so that I still couldn't understand what they are exactly, where I can find them, and if I cannot find them, what I should do so that downstream analyses could proceed (yes, I guess the answer is to add "some" read groups, but what exactly I should add ?)

From what I've found, ID, SM, PL, LB seem important read groups (for GATK at least). But if I am to add these read groups to my BAM files, assuming the files don't have them, can I just assign some dummy names to each of them? Okay, PL probably needs to be specific, like either illumina, solid, or others, but does it matter if I assign them all lowercase or should they be all CAP ?? What about the other RGs ?

For example, if I have only one BAM file to add/replace the read groups, could I simply assign "A", "B", "illumina" and "D" for ID, SM, PL, LB respectively.

And if I have two BAM files, could I simply assign "A1, B1, illumina, D1" for file 1 and "A2, B2, illumina, D2) for file 2?

I found that GATK forum mentioned that dummy info is OKAY, so would A,B,C,D like the examples above be fine ? And what exactly are the purposes for these read groups? If they are so essential, why couldn't they be incorporated by default when running early steps (or even 1st step, e.g. from fastq) of NGS data processing ?

Any input on any of the issues in this question will be greatly appreciated. Thank you.

Answer: Yes, you can assign dummy names for any and all of these. The read group tags are meant to enable grouping of alignment to account for biases due to things like the library preparation, the machine things were sequenced on, etc.

This is mostly useful where you have samples that were each sequenced multiple times, but from different libraries. So then you'd have alignments with the same SM but a different LB. In cases where you just have a single run of each sample, with all samples done in a single batch, then read groups aren't particularly useful.
Question: Hi All,

I was wondering how important is to get rid of exact duplicate Illumina reads before --

1. Before using it for correcting PacBio reads (planning to use ProovRead)
2. Before using it to polish a Pac-Bio only assembly using Pilon (Assembly was done using uncorrected PacBio reads - miniasm)
3. Before using the reads to do a hybrid de-novo-assembly using PBcR

Some of my Illumina libraries have significant amounts of reads duplicated >10 times. What are your recommendations to handle these duplicate reads considering the scenarios mentioned above?

Many thanks in advance!

Answer: It's not a good idea to remove duplicate reads unless your libraries are amplified. If they *are* amplified, and you have reads appearing 10+ times, I highly recommend you change to an unamplified protocol, because you are wasting sequence. And by duplicates, I mean that both read 1 and read 2 of pairs are duplicates... otherwise the pairs are not, in fact, duplicates.

But - if you have a situation in which you are using an amplified library, and duplicate pairs occur, I recommend eliminating all duplicates and replacing them with a single copy of their consensus, in any situation other than quantification (e.g. RNA-seq).
Question: I am new to R but I would like to explore this package. I read a similar posted question about installing this package using a Linux machine, but I was not able to successfully use the information posted there.

The paper describing this tool links to https://github.com/csuzfq/mapsnp_pkg/tree/master/mapsnp, but I am not sure what to do with files from github.

Thank you for any help.

Answer: Next time please post the error when you ask for help and also any material you used.

There's 2 errors your getting and I can replicate them on my system (macos). So the good (or bad) news is that it's not just you. First, the two packages skipped are on Bioconductor, if you google you'll find out how to install them.

Second, I think this is an issue with how devtools is trying to install. Something has gone wrong probably with how the authors have setup the package on github (badly). I suggest you download the source files and install manually like this:

install.packages("~/Downloads/mapsnp_0.1.tar.gz", repos = NULL, type = "source")

This is how I installed the package.

Source files from here: https://github.com/csuzfq/mapsnp_pkg/blob/master/mapsnp_0.1.tar.gz

P.S Make sure you install the two bioconductor packages before installing mapsnp. Then load the two bioc packages first then load mapsnp just to be sure.
Question: Dear Biostars,

I want to prepare the rank file for GSEA analysis based on RNA-seq results that were generated by DESeq2. I have found different recommendations as how to create the pre-ranked gene list. The GSEA site mentions the gene list can be sorted by any value, however, other people have pointed out in this blog that the direction of fold change is important for GSEA analysis. Now, if the genes are sorted only based on their log fold change, a gene with a large fold change but a poor p-value will be ranked higher than a gene with a statistically significant fold change that is smaller in magnitude.

I've also read Mark Zeimann's post about his approach to this issue where he generates a new scoring metric by multiplying the sign of fold change by its inverse p-value: http://genomespot.blogspot.com.au/2014/09/data-analysis-step-8-pathway-analysis.html

He also adds that: "at the top of the list are the genes with "strongest" up-regulation and the bottom of the list are genes with "strongest" down-regulation and genes not changing are in the middle". I am not sure if this is the right assumption for GSEA input file?

I greatly appreciate if you could help me understand this, and explain me your preferred method for creating a GSEA rank file for RNASeq expression results.

Many Thanks,
Noushin

Answer: Hi Noushin, As you note there is no set ranking method for GSEA. There are many alternatives. Some people use fold change. In our group we use signed p-value. It is not a new method, nor did I invent or "generate" it. We have published this extensively and no reviewer has ever commented on the validity. I've included a list of papers in the last 2 years using this method. Now what was your specific question? Cheers, Mark

http://www.ncbi.nlm.nih.gov/pubmed/27157830
http://www.ncbi.nlm.nih.gov/pubmed/27108813
http://www.ncbi.nlm.nih.gov/pubmed/25941940
http://www.ncbi.nlm.nih.gov/pubmed/25491312
http://www.ncbi.nlm.nih.gov/pubmed/24875254
http://www.ncbi.nlm.nih.gov/pubmed/24850797
http://www.ncbi.nlm.nih.gov/pubmed/24732587