cannin/biostars_qa huggingface dataset:
Question: I have paired-end RNA-Seq data - read1 and read2 - stored in the same fastq file.

I'd like to align the reads using tophat.

Do I have to separate the data into two different files before running tophat?

Answer: If you have any pattern in the read name ( like `1:N:####` or `/1` or `_1` etc) you could use the fastq-grep to match the pattern to extract the R1 and R2 into two separate files.

Or a simple Awk patter match will do. something like:

zcat fastq.gz | paste - - - - | awk '{ if $1 ~ /< R1 pattern here>/ { print $1"\n"$2"\n"$3"\n"$4" }' | gzip > Read_1.fastq.gz

zcat fastq.gz | paste - - - - | awk '{ if $1 ~ /< R2 pattern here>/ { print $1"\n"$2"\n"$3"\n"$4" }' | gzip > Read_2.fastq.gz

But they need to be kept in order, if they are not.
Question: Hello,

Our lab would like to use [Homer][1] and install it under the following environment:

- Mac OS X 10.8.5
- Xcode Version 5.0.1 (5A2053)
- Xcode command line tools
- Homebrew
- Python 2.7.2 @ /usr/bin/python

It looks like Homer was successfully installed except for one package, seqlogo. The results of our troubleshooting is listed below. We attempted to contact technical support at Homer and Seqlogo but unfortunately never heard back.

How can we get Homer to see seqlogo? Thanks!

**Troubleshooting:**

```
perl configureHomer.pl -check

Current base directory for HOMER is /Users/Shared/homer/./

Checking for standard utilities and 3rd party software:

Checking for wget... /usr/local/bin/wget
Checking for cut... /usr/bin/cut
Checking for gcc... /usr/bin/gcc
Checking for zip... /usr/bin/zip
Checking for unzip... /usr/bin/unzip
Checking for make... /usr/bin/make
Checking for tar... /usr/bin/tar
Checking for gunzip... /usr/bin/gunzip
Checking for gzip... /usr/bin/gzip
Checking for gs... /usr/local/bin/gs
**Checking for seqlogo..The program seqlogo was not found but is required for making motif logos**
Checking for blat... /usr/local/bin/blat

All auxilary programs found.
```

```
$ seqlogo
Command not found
$ which weblogo
/usr/local/bin/weblogo
weblogo 3.4
```

webLogo installed like so

```
sudo pip install numpy
brew install pdf2svg
sudo pip install weblogo
```

`findMotifs.pl` returns help page (good sign, right?)

```
echo $PATH
/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin:/Users/administrator/Downloads/software/ngsplot/bin:/Users/Shared/homer/bin:/opt/X11/bin
```

**Edit 2014.09.24**

As suggested by Ying W, weblogo 2.8.2 must be used in order for seqlogo to be recognized by HOMER.

[1]: http://homer.salk.edu/homer/introduction/install.html

Answer: Follow the weblogo install directions found here you will need a specific version that might not be the right one installed by pip
Question: Hi,

I've downloaded a set of bams from Illumina's Platinum Genomes experiment (hosted at EBI) and would like to do due diligence and perform a checksum check, but I can't find any published for these files.

Has anyone else downloaded the bam files from ftp://ftp.sra.ebi.ac.uk/vol1/ERA172/ERA172924/bam/ and run `md5sum` on them? These are the CEPH pedigree sequenced to 50x depth.

If so I'd be very grateful for the results to do a comparison (I can publish the checksum values I get once its finished running).

SB

Answer: EBI are actually very good at storing their data and usually have the md5 checksums among the files.

To find the md5 number for your data: look at their data table for the accession number instead of raw FTP access: https://www.ebi.ac.uk/ena/data/view/ERA172924&display=html

Click on "select columns", then on "submitted md5" and voila -- md5 is there. There's a bunch of other metadata there as well.

NB: If you want to parse this from code checkout the [text view][1].

[1]: https://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=ERA172924&result=read_run&fields=study_accession,secondary_study_accession,sample_accession,secondary_sample_accession,experiment_accession,run_accession,tax_id,scientific_name,instrument_model,library_layout,fastq_ftp,fastq_galaxy,submitted_ftp,submitted_galaxy
Question: Hi
I've done PCA for my gene expression data after DEG analysis, I can see that my case and control samples clustered distinctly but **I want to report a p.value for this result**, it worth to notice that I have the coordinate csv file that represents coordinates of samples from different phenotypes (case vs control).

Thanks in advance

Answer: PCA is an exploratory data analysis method. It does not test a null hypothesis and generate a p-value.

If you want to compute a p-value maybe you should try pvclust package in R. It didn't use PCA but a hierarchical clustering and report p-values for each sub-tree

http://stat.sys.i.kyoto-u.ac.jp/prog/pvclust/
Question: I got RNASeq data in several samples. I checked the FastQC, seems the read quality are good (Hiseq 2000). But the problem is many reads are mapped to intronic region, and the regions have no any reference exons there (Refseq, ensembl, gencode). We don't know what they are. We guess the problem happend in library preparation, the concentration was low. Now the data has come out and we can't re-sequencing, so we want to remove the reads mapped to intronic region, is there a method to do that? Or anyone have an idea about the intronic reads. Thanks.

Answer: If you have a bed file of exonic regions, or gtf, something like that, you can use BEDTools to filter your .bam for reads that fall in the desired coordinates, using intersectBed
Question: Hello,
Actually I am doing a tutorial (https://hakyimlab.github.io/psychencode/generate_weights.html) where they indicate that I need a dbSNP150 reference table containing "chromosome, position, ref, alt, rsid, and dbSNPBuildID" information for the hg19 version of the human genome. I am stuck in this step, can anyone help me and indicate where can I find this information?

Answer: Since you are working in R, you can often use annotation that is provided as part of Bioconductor. The latest release of Bioconductor for example contains a copy of dbSNP [version 150][1] and [even newer ones][2].

However, the annotation provided by Bioconductor refers to another reference genome build (hg38 instead of the long outdated version hg19). Older (ancient!) versions of Bioconductor might still contain hg19 coordinates, but better use [LiftOver][3] to convert the latest annotation back to the older reference genome build.

Alternatively, there is also [a hg19 annotation from dbSNP 135][4] in the current build.

[1]: https://bioconductor.org/packages/3.15/data/annotation/html/SNPlocs.Hsapiens.dbSNP150.GRCh38.html
[2]: https://bioconductor.org/packages/3.15/data/annotation/html/SNPlocs.Hsapiens.dbSNP155.GRCh38.html
[3]: https://bioconductor.org/packages/release/workflows/html/liftOver.html
[4]: https://bioconductor.org/packages/3.15/data/annotation/html/FDb.UCSC.snp135common.hg19.html
Question: I'm reading [this excellent paper][1] which describes a methodology to standardize variant benchmarking process. They say that the normal binary classification form (i.e., TP, FP, FN, and statistics derived from these) are not simple for variant calls. So they go on to describe how they do this, in tabular form.

It's also worth noting that this is how they define sensitivity and specificity:

sensitivity (the ability to detect variants that are known to be present or “absence of false negatives”)
and specificity (the ability to correctly identify the absence of variants or “absence of false positives”)

They do not use specificity, opting instead for precision because "precision is often a more useful metric than specificity due to the very large proportion of true negative positions in the genome."

![contingency table][2]

I'm having a hard time understanding why FP is the way it is. I would've thought that all the n/a's in the first column (e,g, "ref/var2" and ref/var3", "var1/var3" etc) would be FP as well.

It's also hard to decipher why the n/a's occur in the rest of the table?
This may have something to do with their comment:

"Note that we have chosen not to include true negatives (or consequently specificity) in our standardized definitions. This is due to the challenge in defining the number of true negatives, particularly for indels or around complex variants."

Are the n/a's representing true negatives? It is all so confusing.

[1]: https://www.biorxiv.org/content/biorxiv/early/2018/05/24/270157.full.pdf
[2]: https://www.biorxiv.org/content/biorxiv/early/2018/02/23/270157/T1.medium.gif

Answer: > It's also hard to decipher why the n/a's occur in the rest of the table?

I think it is because those scenarios can't happen. If you compar two individuals, and the Query has var1/var3 then there has to be at least a var2 in your comparison, otherwise you can't have var3.

If Query = GT:1/3 then the only scenarios for Truth are GT:1:2, in a vcf file.
Question: Hi all.

I have detected somatic mutation and would like to know its significance for clinical interpretation.

I used MutationMapper in cBioPortal and then got a set of annotated mutations with plot as shown in below.

As you can see, there are a protein domains in the plot with yellow color and all the mutations are on its domain.

and my questions is below:

The segments that do not belong to yellow domains whose color is grey have no meaning for clinically annotating mutation? so I discard it unconsciously?

Answer: Your question is an important one without a definitive answer. In terms of clinical annotation, protein domain is just one of MANY, MANY potential features for prioritization. In practice, I would not discard mutations that are not in a known domain.
Question: I've got a list of 237 Ensembl protein IDs (e.g. ENSP00000493027), and I'm trying to convert them to UniProt accession numbers so that I can retrieve their REST text entry (e.g. https://rest.uniprot.org/uniprotkb/A0A286YF28.txt).

Officially, UniProt says to do it [this way][1]:

import urllib.parse
import urllib.request

url = 'https://www.uniprot.org/uploadlists/'

params = {
'from': 'ACC+ID',
'to': 'ENSEMBL_ID',
'format': 'tab',
'query': 'P40925 P40926 O43175 Q9UM73 P97793'
}

data = urllib.parse.urlencode(params)
data = data.encode('utf-8')
req = urllib.request.Request(url, data)
with urllib.request.urlopen(req) as f:
response = f.read()
print(response.decode('utf-8'))

However, when I run their example code verbatin, I get:

urllib.error.HTTPError: HTTP Error 405: Not Allowed

Along with some tracebacks. What's going on? It happens when I pass my own data too. Never used UniProt programmatically before.

Thanks in advance for any help anyone can give.

[1]: https://www.uniprot.org/help/api_idmapping

Answer: Here is an example how to access UniProt's REST with Python 3 with the `requests` package (`pip install requests`).

import json
import requests
import time

URL = 'https://rest.uniprot.org/idmapping'
IDS = ['P40925', 'P40926', 'O43175', 'Q9UM73', 'P97793']

params = {
'from': 'UniProtKB_AC-ID',
'to': 'Ensembl_Protein',
'ids': ' '.join(IDS)
}

response = requests.post(f'{URL}/run', params)
job_id = response.json()['jobId']
job_status = requests.get(f'{URL}/status/{job_id}')
d = job_status.json()

# Make three attemps to get the results
for i in range(3):
if d.get("job_status") == 'FINISHED' or d.get('results'):
job_results = requests.get(f'{URL}/results/{job_id}')
results = job_results.json()
for obj in results['results']:
print(f'{obj["from"]}\t{obj["to"]}')
break
time.sleep(1)

Output:

P40925	ENSP00000233114.8
P40925	ENSP00000410073.2
P40925	ENSP00000446395.2
P40926	ENSP00000327070.5
P40926	ENSP00000408649.2
O43175	ENSP00000493175.1
O43175	ENSP00000493382.1
Q9UM73	ENSP00000373700.3
P97793	ENSMUSP00000083840
Question: My experiment set up: 2 samples from WT mice, 2 samples from KO mice, all sequenced with 10x 3' scRNA-seq. The cell populations sorted for sequencing are the same.

Seurat offers the anchor transfer method to perform dimension reduction and clustering, but differential expression is still performed upon normalized, untransformed data matrices.

Therefore, I was wondering if there are methods available to perform batch correction that can be used in differential expression. I would appreciate it if you could share your experience in dealing with this type of situation.

Thanks a lot~

Answer: I personally like to aggregate cells per genotype and cluster into pseudobulks and then simply include the batch information into the design as we do in any normal RNA-seq setup such as `~batch+pseudobulkCluster`. Since you have replicates per genotype this comes then down to a normal 2 vs 2 comparison. That having said the few samples I analyzed were always made in a way that there was little batch effect so we processed one replicate of each genotype on the same day without any inter-platform or inter-species comparisons which probably is beyond simply including batch as a covariate. But for me that strategy worked out well so far. I personally use `sumCountsAcrossCells` from the scater package (assumes the SingleCellExperiment format) to get the pseudobulks and then edgeR or DESeq2 for the DE.
Question: I have a dataset with a missing column header for first column. How to name it so that I can find it using

a$
![screenshot][1]
[1]: https://image.ibb.co/nNv7Xn/column_name.jpg

I used following command but got error

a$miRNAs<-row.names(a)

Warning message:
In a$miRNAs <- row.names(a) : Coercing LHS to a list

How to find the solution?

Answer: It is rownames not a column. Since you have mentioned dplyr, we can convert it to columns using packages from tidyverse, here is an example:

# example dataset with rownames
a <- mtcars[1:3, 1:3]
a
# mpg cyl disp
# Mazda RX4 21.0 6 160
# Mazda RX4 Wag 21.0 6 160
# Datsun 710 22.8 4 108

library(dplyr)
library(tibble)

a <- a %>%
rownames_to_column(var = "myName")

a
# myName mpg cyl disp
# 1 Mazda RX4 21.0 6 160
# 2 Mazda RX4 Wag 21.0 6 160
# 3 Datsun 710 22.8 4 108
Question: It might be a very stupid question for many of you but, since it's my first variant calling, I didn't figure it out yet.

I have **mpileup**'ped two bam files from two samples, then I filtered the results with **vcfutils.pl** and called the genotypes with **bcftools** **call**. Now I have a **VCF** file containing what I want, but I managed to analyze differences only between the two samples and the reference (which is an assembly coming from a line that is *different* *from both* *samples*).

I have the variants between my two samples and the assembly, but what if I want to detect the differences **between** the two samples? I did it with awk / sed / cut and other command line tools but is there maybe a better and more straight-forward way to do that? Perhaps using bcftools? Up to know, I didn't find it.

Any suggestion appreaciated!

EDIT:
I did it with bcftools gtcheck as well, that works fine. I am asking if there are **other** reliable tools to test!

EDIT 2:
As this post has many views now, I guess it will be useful for everyone to know that I used **bcftools isec** and it worked brilliantly.

Answer: depending on what you're looking for, like simple counts for instance, awk/perl solutions could be enough. if you need a deeper description, `bcftools stats` would give you some interesting details.

if your 2 samples are in a single vcf file I would suggest to split them first

bcftools view -Oz -c1 -s sample1 joined.vcf.gz > sample1.vcf.gz
bcftools view -Oz -c1 -s sample2 joined.vcf.gz > sample2.vcf.gz

and then compare them

bcftools stats sample1.vcf.gz sample2.vcf.gz > joined.stats.txt
Question: **Can I transfer this first format to the second one just by basic shell procession or awk or sed on linux? This is a toy example:**

This kind of text file is what I have, three cols, col2 and col3 like range, left close and right open,
chr1 0 2 0
chr1 2 6 1.5
chr2 0 3 0
chr2 3 10 2.1

Transfer to describe each position as:
chr1 0 0
chr1 1 0
chr1 2 1.5
chr1 3 1.5
chr1 4 1.5
chr1 5 1.5
chr2 0 0
chr2 1 0
chr2 2 0
chr2 3 2.1
...
chr2 9 2.1

Someone has idea how to solve? Thanks!!!

Answer: awk '{B=int($2);E=int($3);for(i=B;i
Question: I need to calculate the median in the fastq file, how can I do this with awk or bioawk?

Answer: Straight up `awk` soliution would be:

cat *.fq | awk 'NR % 4 == 0 { print length($1) } ' | awk ' { a[i++]=$1; } END { print a[int(i/2)]; }'

I would do it with `bioawk` and `datamash`

cat *.fq | bioawk -c fastx '{ print(length($seq)) }' | datamash median 1

you can also implement a better median in awk as shown here:

* https://stackoverflow.com/questions/6166375/median-of-column-with-awk

that works like this:

cat *.fq | bioawk -c fastx '{ print(length($seq)) }' | bioawk ' { a[i++]=$1; } \
END { x=int((i+1)/2); if (x < (i+1)/2) print (a[x-1]+a[x])/2; else print a[x-1]; }'
Question: Hi Everyone,

We have commissioned RNA-seq and analysis by a company, which provided us with raw fastq files, BAM files, and a count matrix. They used hard clipping and Tophat for the alignment to GRCm38/Mm10. I have attempted to recreate their analysis with HISAT2 (same reference genome), using simply the default parameters and no separate trimming/clipping. I have used samtools to convert the SAM files to BAM files and compared the results from the company's analysis ("Tophat") with my own ("HISAT2") using IGV. The results are very confusing to me. The majority of genes I have (randomly) inspected look highly similar between both sets of BAM files. See this example gene (Tophat in blue, HISAT2 in red):

[Tnf][1]

So far, so good. However, there are also multiple instance where one analysis picked up good reads, while the other did not. This is true in both directions. See these two example genes:

[Gapdh][2]

[Il12b][3]

And, finally, there are some genes in which one alignment just looks weirdly skewed. For instance:

[Ubc][4]

Does anyone know what might account for these differences? Or which alignment I should use for downstream analysis? I'd be grateful for any feedback!

Thomas

[1]: https://ibb.co/ksUQsS
[2]: https://ibb.co/b1emz7
[3]: https://ibb.co/jj9QsS
[4]: https://ibb.co/nO9mz7

Answer: Apparently, this problem is inherent to IGV. When we ran the same files that showed no reads on a different operating system, everything looked just fine. There are, of course, still differences between the Tophat and the HISAT2 alignments, but nothing is missing altogether. I will try to switch to the UCSC Genome Browser, as suggested by WouterDeCoster.
Thank you all for your input!
Question: Hello All,

I have read counts from RNA seq data in row and columns. I want to quantile normalized them in R. I have following code. This gives me the normalized values. However, the output is a matrix. I want the output with row name and column name so that I can perform PCA on it.

data <- read.csv("data.csv",header=T)
head(data)
data_mat <- as.matrix(data[,-1])
head(data_mat)
data_norm <- normalize.quantiles(data_mat, copy = TRUE)

Could someone help me to get that? Thank you in advance.

Answer: Try this (note the extra line; also use `data.matrix`, not `as.matrix`):

data <- read.csv("data.csv",header=T)
head(data)
rownames(data) <- data[,1]
data_mat <- data.matrix(data[,-1])
head(data_mat)
data_norm <- normalize.quantiles(data_mat, copy = TRUE)
Question: So I'm trying to understand the --shift and --extsize parameters in MACS2. I have inherited an ATAC-seq that used the following MACS2 command.

macs2 callpeak -t ATAC_sample-1.bam --nomodel --shift -100 --extsize 200 -g 1.5e9 -f BAM

I'm trying to see how I can apply these principles to a Rscript I have that converts my raw BAM file to a coverage file? The full script can be looked at [here](http://rpubs.com/achitsaz/98857) but the main lines I'm looking at are:

aln <- as(aln, "GRanges") # Converts each read mapping to GRanges coord

aln <- resize(aln, 150) # Extended the reads to the fragment length (4 previous exp)

cov <- coverage(aln) # Get Coverages nucleotide

The definitions for extsize and shift can be seen in the [documentation](https://github.com/taoliu/MACS#--extsize) but I'll copy it for reference.

> --extsize

> >While '--nomodel' is set, MACS uses this parameter to extend reads in 5'->3' direction to fix-sized fragments. For example, if the size of binding region for your transcription factor is 200 bp, and you want to bypass the model building by MACS, this parameter can be set as 200. This option is only valid when --nomodel is set or when MACS fails to build model and --fix-bimodal is on.

>--shift

>>Note, this is NOT the legacy --shiftsize option which is replaced by --extsize! You can set an arbitrary shift in bp here. Please Use discretion while setting it other than default value (0). When --nomodel is set, MACS will use this value to move cutting ends (5') then apply --extsize from 5' to 3' direction to extend them to fragments. When this value is negative, ends will be moved toward 3'->5' direction, otherwise 5'->3' direction. Recommended to keep it as default 0 for ChIP-Seq datasets, or -1 * half of EXTSIZE together with --extsize option for detecting enriched cutting loci such as certain DNAseI-Seq datasets. Note, you can't set values other than 0 if format is BAMPE or BEDPE for paired-end data. Default is 0.

>Here are some examples for combining --shift and --extsize:

>>EXAMPLE 1: To find enriched cutting sites such as some DNAse-Seq datasets. In this case, all 5' ends of sequenced reads should be extended in both direction to smooth the pileup signals. If the wanted smoothing window is 200bps, then use '--nomodel --shift -100 --extsize 200'.

So to summarize my question is how should I resize my raw reads so that the coverage will reflect the --shift and --extsize used by MACS2? From reading the documentation it seems I should shift the reads 100 bp in both directions because it is extending the reads 200 bp 5' -> 3' (--extsze) and back 100 bp with (--shift).

Answer: You'll need to `endoapply()` a function that does the `shift()` according to the strand, since that's otherwise ignored. Then `resize(aln, 200)`. Alternatively, use `idx = which(strand(aln) == '+')` for each strand and directly adjust the start accordingly (then `resize()`).
Question: Hello,

I have microarray data from the chip [HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version], and I am trying to annotate the probe IDs in R using the hugene10sttranscriptcluster.db package and the following function:

annotatedTopTable <- function(topTab, anotPackage)
{
topTab <- cbind(PROBEID=rownames(topTab), topTab)
myProbes <- rownames(topTab)

thePackage <- eval(parse(text = anotPackage))

geneAnots <- AnnotationDbi::select(thePackage, myProbes, c("SYMBOL", "ENTREZID", "GENENAME"))

annotatedTopTab<- merge(x=geneAnots, y=topTab, by.x="PROBEID", by.y="PROBEID")
return(annotatedTopTab)
}

topAnnotated_Condition<- annotatedTopTable(topTab_Condition,
anotPackage="hugene10sttranscriptcluster.db")

topAnnotated_Condition

Is hugene10sttranscriptcluster.db the right package to use?

Thanks

Answer: Yes, for most use-cases, that package is correct. It will depend on how the original array data was processed.

You only need to do this, by the way:

require(hugene10sttranscriptcluster.db)
mapIds(
hugene10sttranscriptcluster.db,
keys = probes,
column = 'SYMBOL',
keytype = 'PROBEID')

select(
hugene10sttranscriptcluster.db,
keys = probes,
column = c('SYMBOL', 'ENTREZID', 'ENSEMBL'),
keytype = 'PROBEID')

Kevin
Question: Hi folks:

I am trying to get raw RNASeq data from ENA and wonder how to know the length of the sequences in those FASTQ files without downloading them?

In the query design I follow their [guide][1], and use the "`nominal_length>XXX"` in the query, but this filtering mechanism doesn't work...

Thanks!

[1]: http://www.ebi.ac.uk/ena/browse/search-rest

Answer: Get the count of total bases and number of reads from ENA, and do below calculation,

Read length = Total Bases/ number of reads