cannin/biostars_qa huggingface dataset:
Question: Hi all,

I am working with RNA-seq data of different drug treated human cell lines sequenced on Illumina (2X150 bp chemistry). Data generated is 11-14 Gb.

I am using `HISAT2` for alignment on human genome hg38 build downloaded from [Ensembl database][1] using default parameters. To my surprise, the alignment percentage is very low (`~2 %`) for all the samples.

I have the following observations regarding the data quality :

1. Data quality (phred score) is good i.e. in the range of `32-40`!
2. Illumina adapters were already removed using trimmomatic. No other contaminants.
3. The fastqc metric '[Per base sequence content][2]' fails for all samples. More specifically, A,T,G,C % is deviating considerably at the first 10 base positions.

**Questions**

1. What could possibly be wrong? Any other QC stuff I need to re-consider?
2. More importantly, this is cell-line data, so shall I consider mapping on specific cell line reference OR is it okay if I have just the hg38 standard human reference?

[1]: ftp://ftp.ensembl.org/pub/release-90/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz
[2]: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/4%20Per%20Base%20Sequence%20Content.html

Answer: I've seen cell lines being contaminated or swapped by other lines - in the "good" cases by cell lines from other species (good because when this happens it is really easy to detect). Your low mapping rate suggests you are mapping to the wrong species, [Centrifuge][1] is really handy to detect this kind of contamination - the problem is to download or build a nt index.

[1]: http://www.ccb.jhu.edu/software/centrifuge/
Question: I've been using fgsea to run some enrichment analysis and one of the outputs from `fgsea()` is *nMoreExtreme*.

I've seen from the vignette that the *nMoreExtreme* is:

> a number of times a random gene set had a more extreme enrichment score value

I was wondering what this actually means?

Like what is a random gene set? Is it just a random selection of genes? But then would nMoreExtreme be different every time fgsea is run?

Is it another measure of how reliable the enrichment analysis is?

I am guessing that a value of 0 or 1 is good? And a value of a few thousand is bad?

Could I use it alongside a padj value to determine if the results are significant or is it not as important as padj?

Thanks in advance.

Answer: GSEA P-values are empirically calculated by sampling random gene sets. `nMoreExtreme` is indeed a number of times a random gene set had a more extreme enrichment score value and is used in nominal GSEA P-value calculation, which is equal `nMoreExtreme` normalized by number of random gene sets that have the same sign of enrichment score. As it's random, it can have different values from run to run. If it's zero it means that the calculated P-value is bounded by the number of permutations, not the real P-value, and the result can be huge overestimation (or not).

These said, in the recent release the algorithm has be changed to fgsea-multilevel, now there are no more `nMoreExtreme` values and the estimation errors are reported explicitly.
Question: Is there a tool or maybe a method in R to produce a line graph of the coverage of mapped resequencing, for each chromosome of the mapped bam file? I'm looking to compare the coverage across the genome of mapped reads using two different mapping software.

Answer: Do a simple loop with bash & samtools:

cut -d ' ' -f 1 reference.fa.fai |\
while read CHROM
do
for BAM in *.bam
do
echo -n -e "${CHROM}\t${BAM}\t"
samtools view -F 4 -c ${BAM} ${CHROM}
done
done
Question: For each of several genomes, apart from the already available fasta sequence and associated GFF3 annotation files, I have also generated 5 additional GFF files for start-stop coordinates of 5 additional types of genomic features.

My goal is four-fold in the context of these 6 GFF files and 1 genomic DNA sequence.

1. **Explore** where two or more of these genomic features overlap / intersect / co-localize - I am doing this via text manipulation, using bedtools overlap or bedtools intersect. But I want to visualize 2 or more tracks , and just text-based calculation of intersection is not satisfying during exploration phase. I want to visualize it.

2. Generate **high-detail images** (with flexibility of color / shapes like that in IGB or gff2ps) for **small intervals**, looking at specific and most interesting cases.

3. Generate **overview images** for **entire** chromosomes or even a **genome**, for overlap across these 6 types of genomic features , without confusing or overwhelming the reader.

4. Finally, to request advice on which tool and/or which **statistical test** to perform for verifying whether the observed physical **co-localization** of any 2 of the 6 types of genomic features is random or **non-random**. And if latter is true, are there more sophisticated tests to examine physical distribution of genomic loci types, relative to one another? And are there tests that can examine more than 2 types of genomic features at a time?

The rather old thread at https://www.biostars.org/p/363/ discusses answers to questions 1 - 3 above, but I am curious to know if there are better / updates tools for my goals than the ones I mentioned above or at the link (bedtools or bedOps, IGB, GBrowse, GFF2PS). Thanks!

Answer: You could do this (I think) with [Artemis][1] and/or the [Artemis Comparison Tool][2].

Load up the sequence, and read in the different annotation files and you can view them all together in the same window I believe.

It probably doesn't solve *all* of your requests though, but worth a look I think.

[1]: https://www.sanger.ac.uk/science/tools/artemis
[2]: https://www.sanger.ac.uk/science/tools/artemis-comparison-tool-act
Question: Hi,
I'll try to explain my problem the best I can.
I have a script that runs multiple bowtie2 alignments. The problem is that it continues running the next alignment even if the current one fails to run.

**script.sh**

#!/bin/bash

bowtie2 --threads 32 -x sp_idx -U reads1.fastq | samtools view -Sbu - | samtools sort -@16 -m 4G > 1_sorted.bam
bowtie2 --threads 32 -x sp_idx -U reads2.fastq | samtools view -Sbu - | samtools sort -@16 -m 4G > 2_sorted.bam
bowtie2 --threads 32 -x sp_idx -U reads3.fastq | samtools view -Sbu - | samtools sort -@16 -m 4G > 3_sorted.bam
bowtie2 --threads 32 -x sp_idx -U reads4.fastq | samtools view -Sbu - | samtools sort -@16 -m 4G > 4_sorted.bam

What I want is that if there's an error while running Bowtie2, the script stops running and not keep running to the end. What's happening is that if line 2 fails, the script continues with the alignment of sample 3 and then 4. So why doesn't it stop running if line 2 is errors out?

Thanks!

Answer: That's how shell scripts operate. If you expect that your shell script should stop as soon as an error (exit code != 0) is encountered, use `set -eo pipefail` right at the top after the interpreter (`!#/bin/bash`) line.

#!/bin/bash
set -eo pipefail

bowtie2 --threads 32 -x sp_idx -U reads1.fastq | samtools view -Sbu - | samtools sort -@16 -m 4G > 1_sorted.bam
bowtie2 --threads 32 -x sp_idx -U reads2.fastq | samtools view -Sbu - | samtools sort -@16 -m 4G > 2_sorted.bam
bowtie2 --threads 32 -x sp_idx -U reads3.fastq | samtools view -Sbu - | samtools sort -@16 -m 4G > 3_sorted.bam
bowtie2 --threads 32 -x sp_idx -U reads4.fastq | samtools view -Sbu - | samtools sort -@16 -m 4G > 4_sorted.bam

More reading:

- The `set` builtin command: https://www.gnu.org/software/bash/manual/html_node/The-Set-Builtin.html (check out the `-e` [which is what you're looking for] and the `-o pipefail` [which is a special case of your requirement where a command that is part of a pipe fails])
Question: Hello,

I need your help. I want to compare two dataframes with genomic features like this:

More precisely, frame 2 compared to frame 1

chr start end name chr start end name
1 2 5 A 1 1 3 AA
2 10 15 B 2 9 16 BB
3 27 30 C 3 28 29 CC

As a result I want this:

1. Completely overlapping:

```
chr start end name chr start end name
2 10 15 B 2 9 16 BB
```

2. Partial overlap:

```
chr start end name chr start end name
1 2 5 A 1 1 3 AA
```

3. and within the range:

```
chr start end name chr start end name
3 27 30 C 3 28 29 CC
```

Are the R packages `IRanges` and `GenomicRanges` suitable for such analysis? Or do I have to write some > < commands?

Answer: Using the GenomicRanges library in R:

library(GenomicRanges)
x1 = read.table(text="chr start end name
1 2 5 A
2 10 15 B
3 27 30 C",header=T)

x2 = read.table(text="chr start end name
1 1 1 3 AA
2 2 9 16 BB
3 3 28 29 CC",header=T)

# Make GRanges object
gr1 = with(x1, GRanges(chr, IRanges(start = start, end = end, names = name)))
gr2 = with(x2, GRanges(chr, IRanges(start = start, end = end, names = name)))

# Completely overlapping
type1 = findOverlaps(query = gr1, subject = gr2, type = "within")
type1.df = data.frame(x1[queryHits(type1),], x2[subjectHits(type1),])
type1.df
chr start end name chr.1 start.1 end.1 name.1
2 2 10 15 B 2 9 16 BB

# Within range
type3 = findOverlaps(query = gr2, subject = gr1, type = "within")
type3.df = data.frame(x1[subjectHits(type3),], x2[queryHits(type3),])
type3.df
chr start end name chr.1 start.1 end.1 name.1
3 3 27 30 C 3 28 29 CC

# Partial Overlaps only (no complete overlaps or within range overlaps)
type2 = findOverlaps(query = gr1, subject = gr2, type = 'any')
type2.df = data.frame(x1[queryHits(type2),], x2[subjectHits(type2),])
x = rbind(type1.df, type2.df, type3.df)
type2.df = x[!(duplicated(x) | duplicated(x, fromLast = TRUE)), ]
type2.df
chr start end name chr.1 start.1 end.1 name.1
21 1 2 5 A 1 1 3 AA

You will get three different data frames as per your question. Based on your example, I am assuming type 1 and type 3 are opposites.
Question: I'm running the dropseq pipeline and there is a part where the samfile gets sorted. It looks like there is an option to either sort by coordinates or by queryname. Is there a benefit to either of these?

Answer: Reading a file sequentially is faster than random access, and keeping in memory just the information necessary for some calculation is more efficient than keeping the whole file. Some tasks are more easily performed depending on how the bam is sorted, because the bam can be read sequentially and just part of the data need to be kept in memory.

For example, marking duplicates (which, for paired reads, is done by looking at 5' mapping positions of both reads) is a lot easier for bams sorted by position, because you guarantee the reads physically closer inside the bam are also close on the genome. If they weren't, one would need to scan the whole file to build a hash of reads per position in order to mark duplicates.

Conversely, counting reads mapping to features it easier for name-sorted files, as read pairs are next to each other, and secondary / supplementary alignments are next to primary alignments. Again, if they weren't, one would need to scan the whole file to build a hash of reads names per feature mapped.

Of course, most immediately as an end-user, one has to pay attention to which settings are necessary and which sorting order is expected by the tool of choice.
Question: Hi,
I tried to follow this [instruction][1] on how to retrain Augustus with BUSCO. I ran BUSCO:

run_BUSCO.py -i NbV1Ch07.fasta -o augustus_busco -l /assembly/busco-embryophyta_odb9/ -m geno --cpu 1 --long -sp coyote_tobacco --augustus_parameters='--progress=true'

and it created a folder `/augustus/run_augustus_busco/augustus_output/retraining_parameters` with the following content:

BUSCO_augustus_busco_3750827137_weightmatrix.txt
BUSCO_augustus_busco_3750827137_parameters.cfg.orig1
BUSCO_augustus_busco_3750827137_parameters.cfg
BUSCO_augustus_busco_3750827137_metapars.utr.cfg
BUSCO_augustus_busco_3750827137_metapars.cgp.cfg
BUSCO_augustus_busco_3750827137_metapars.cfg
BUSCO_augustus_busco_3750827137_intron_probs.pbl
BUSCO_augustus_busco_3750827137_igenic_probs.pbl
BUSCO_augustus_busco_3750827137_exon_probs.pbl

Which file should I use to run Augustus trained with BUSCO and how the command looks like?

Thank you in advance,

[1]: https://www.biostars.org/p/261203/#295503

Answer: You copy the files from `/augustus/run_augustus_busco/augustus_output/retraining_parameters` to your `augustus` installation directory under `/augustus-install-directory/config/species/BUSCO_augustus_busco_3750827137/`

usage:
augustus [parameters] --species=BUSCO_augustus_busco_3750827137 queryfilename

or if you want to output to be redirected to a file (see also outfile parameter):
augustus [parameters] --species=BUSCO_augustus_busco_3750827137 queryfilename > output.gff
Question: Hi, all!

I'm trying to download peptide sequence using `biomaRt::getSequence()` function.

The code I used is like below,

seq = getSequence(id = c("ENSP00000367848"),
type = "ensembl",
seqType = "peptide",
mart = mart)

When I ran it, I got a following error.

Error in getSequence(id = c("ENSP00000367848"), type = "ensembl", seqType = "peptide", :
Invalid type argument. Use the listFilters function to select a valid type argument.

According to its vignette, `"Supported types are hugo, ensembl, embl, entrezgene, refseq, ensemblTrans and unigene."`

Why am I getting this error? Can you guys help me?

Thank you!

Answer: The `type` you need is `emsembl_peptide_id`. The list of supported types are not valid as exact string values, I think.
Question: I want to programmatically list all the aliases for a whole number of gene names. I am using Entrez/NCBI and the "gene" database to do this.

For eg:

When I search the ncbi gene database for XRCC1, I get the results

http://www.ncbi.nlm.nih.gov/gene/?term=XRCC1

In the browser NCBI gives you a nice table with the "Name", "Geneid" , "Location" and importantly "Aliases". How do I get these aliases from a Biopython or Entrez eutils query programmatically?

Currently I can get a bunch of information from a function like:

```
from Bio import Entrez

def main(infile_name):
Entrez.email = "myemail@myworkplace.com"
# How do I setup this query to get Aliases?
query = Entrez.esearch(db="gene",term=input_term,retmode="text")
# Process this result and extract Aliases?
result = Entrez.read(query)
print result
```

Thanks for your help

Answer: You might try using the relatively new [Entrez Direct tools][1] that can be run from the command-line. Here is an example command for pulling the aliases for BRCA1[human].

```
esearch -db gene -query "(BRCA1[Preferred Symbol]) AND 9606[Taxonomy ID]" | \
efetch -format xml | \
xtract -pattern Entrezgene \
-element Gene-ref_locus \
-sep "|" -element Gene-ref_syn_E
```

This results in bar-separated aliases following the official symbol.

BRCA1 IRIS|PSCP|BRCAI|BRCC1|FANCS|PNCA4|RNF53|BROVCA1|PPP1R53

[1]: http://www.ncbi.nlm.nih.gov/books/NBK179288/
Question: Dear BioStars Leaders,

I was under the impression that the nucleotide sequence of an mRNA from Refseq would match identically to the Reference Genome assembly. Please correct me if my understanding is incorrect.

At this location : chr1:877831 (GRCh37), the reference nucleotide is "T"

https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr1%3A877831%2D877831&hgsid=581652983_N7DHUPIeX4dCBoDQHXu8dfrKaa0U

The transcript(gene) at this above genomic location is NM_152486.2(SAMD11).

The genomic coordinate chr1:877831 corresponds to following positions in transcript NM_152486.2:

1. Transcript Position : 1107 (starts from the first utr)
2. Codon Position(c) : 1027 (starts from first start codon)

I looked up for the nucleotide at position 1107 in refseq(NM_152486.2), and it is "C".

https://www.ncbi.nlm.nih.gov/nuccore/NM_152486.2?report=gbwithparts&log$=seqview

I have plugged in a mutation "chr1:g.877831T>C" in Mutalyzer(Position Convertor tool) and it translates it to "NM_152486.2:c.1027T>C"
https://www.mutalyzer.nl/position-converter?assembly_name_or_alias=GRCh37&description=chr1%3Ag.877831T%3EC

After that, I plugged in the converted transcript position "NM_152486.2:c.1027T>C" in Mutalyzer(Name Checker tool), and then is gives an error "T not found at position 1107, found C instead".
https://www.mutalyzer.nl/name-checker?description=NM_152486.2%3Ac.1027T%3EC

I am really confused how the reference genome says the nucleotide is "T" , but the mRNA says it is "C". It would be great if someone could explain this?

Addition to my original question : There is a SNP record for this mutation in dbSNP : rs6672356
https://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=6672356

Thanks,
gsr

Answer: RefSeq can be different from the reference genome. These two have a completely different origin. RefSeq is a curated database of transcripts, while reference genomes are results of a genome assembly. Sets of reads/methods are very different. RefSeq is a stable thing, It does not depend on the reference update and this is why HGVS notation for mutation reporting tells to use RefSeq instead of reference genome.

In my practice, I have seen many differences between RefSeq and genomic references. The craziest and most unexpected difference was the difference is single nucleotide delition in the coding region!
Question: I have one large text table, tab separated. The first row is a header. I then have a second text file that contains a subset of the headers in the first file. I want to extract all the columns of the first file whose header is contained in the list given in the second file. Here an example of the inputs and the desired output:

`DATA.TXT`

```
head0 head1 head2 head3 head4
1 25 1364 22 13
2 10 215 1 22
```

`LIST.TXT`

```
head0
head4
```

Desired output:

```
head0 head4
1 13
2 22
```

Answer: Just with a simple google search a hundred of questions/answers related to this topic can be easily found, specially, in **stackoverflow**: http://unix.stackexchange.com/questions/114829/display-column-that-has-a-specific-string

I'd suggest you to try to solve this kind of problems searching on internet and trying yourself before asking. In this case, you can get what you want with a little modification to the solution proposed on the link:

awk -F'\t' 'NR==FNR{arr[$1]++;next}{for(i=1; i<=NF; i++) if ($i in arr){a[i]++;}} { for (i in a) printf "%s\t", $i; printf "\n"}' file2 file1
Question: I'm trying to use `Python` and `Biopython` to fetch metadata for a given assembly identifier. In this case I'm looking for GCF_000005845.2

from Bio import Entrez

# GCF_000005845.2
id_ecoli = "GCF_000005845.2"
esummary_handle = Entrez.esummary(db="assembly", id=id_ecoli, report="full")
record = Entrez.read(esummary_handle, validate=False)
record
# DictElement({'DocumentSummarySet': DictElement({'DocumentSummary': []}, attributes={'status': 'OK'})}, attributes={})

This is the type of data I'm looking for below:
https://www.ncbi.nlm.nih.gov/assembly/GCF_000005845.2/
![enter image description here][1]

I could make a HTML scraper but I don't want to reinvent the wheel if there is already something available.

[1]: https://i.imgur.com/9yr2L6C.png

Answer: If JSON/XML output will be useful to you, the following script can be used.

#!/usr/bin/python

from Bio import Entrez
import json

#Increase query limit to 10/s & get warnings
Entrez.email = ""
#Get one from https://www.ncbi.nlm.nih.gov/account/settings/ page
Entrez.api_key=""

term="GCF_000005845.2"
#Finds the ids associated with the assembly
def get_ids(term):
ids = []
handle = Entrez.esearch(db="assembly", term=term)
record = Entrez.read(handle)
ids.append(record["IdList"])
return ids

#Fetch raw output
def get_raw_assembly_summary(id):
handle = Entrez.esummary(db="assembly",id=id,report="full")
record = Entrez.read(handle)
#Return individual fields
#XML output: https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=assembly&id=79781&report=%22full%22
#return(record['DocumentSummarySet']['DocumentSummary'][0]['AssemblyName']) #This will return the Assembly name
return(record)

#JSON formatted output
def get_assembly_summary_json(id):
handle = Entrez.esummary(db="assembly",id=id,report="full")
record = Entrez.read(handle)
#Convert raw output to json
return(json.dumps(record, sort_keys=True,indent=4, separators=(',', ': ')))

#Test
for id in get_ids(term):
#print(get_raw_assembly_summary(id)) #For raw output
print(get_assembly_summary_json(id)) #JSON Formatted
Question: Dear colleagues,

Does somebody know, if there are specialized databases for Hi-C raw data (paired-end reads)? Of course, there are many raw Hi-C datasets on GEO, but I failed to find any specialized Hi-C database on the Web.

Answer: I wouldn't expect anyone to duplicate storage of raw reads beyond what SRA and GEO do. But as far as aggregating multiple experiments, the [WashU Epigenome Browser][1] has long-range chromatin interaction experiment tracks available. Currently there are 73 tracks available for human, 27 each for mouse and fly.

[1]: http://epigenomegateway.wustl.edu/browser/
Question: I have a FASTA file of gene's like this

>CS1RZP_12510_67 # 64357 # 65316 # -1 # ID=1311_67;partial=00;start_type=ATG;rbs_motif=3Base/5BMM;rbs_spacer=13-15bp;gc_cont=0.675
ATGGCCGTCGTCACCATGAAGCAGATGCTTGATTCCGGCGTGCATTTCGG
GCATCAGACCCGCCGCTGGAACCCGAAGATGAAGCGCTACATCCTGACCG
AGCGCAACGGGATCTACATCATCGACCTGCGGCAGACGCTCACCTACGTC

I would like to extract the corresponding lines for these specific genes from a GFF3 file which contains info on other genes as well which aren't in the above FASTA.

CS1RZP_12510	Prodigal_v2.6.3	CDS	3	617	134.5	-	0	ID=1311_67;partial=00;start_type=ATG;rbs_motif=3Base/5BMM;rbs_spacer=13-15bp;gc_cont=0.675;conf=100.00;score=133.83;cscore=137.59;sscore=-3.75;rscore=-3.31;uscore=-4.34;tscore=4.54;
CS1RZPR4D_68	Prodigal_v2.6.3	CDS	1000	1218	56.4	+	0	ID=1_2;partial=01;start_type=ATG;rbs_motif=AGGA;rbs_spacer=5-10bp;gc_cont=0.612;conf=100.00;score=56.40;cscore=44.46;sscore=11.94;rscore=7.30;uscore=0.10;tscore=4.54;

I am looking for an AWK solution. Grep using -wf is extremely slow

Answer: Figured it out with a little help from a friend, if anyone finds this useful

awk '/^>/,NR==FNR{split(substr($1,2),a,"_");split($9,b,";");HEADER[a[1]"_"a[2]]=b[1]} \
{split($9,c,";");if($1 in HEADER && c[1] == HEADER[$1]) print $0}' sample.fa sample.gff
Question: Hello

I performed a ChIP Seq experiment and performed these steps MACS, CEAS, Gene Ontology Enrichment with GREAT .

I have some enriched terms as well as enriched pathways such as PTEN Signalling pathway, WNT pathway etc.

From my information how can I find the genes which are up or down regulated and the effects they are having on these pathways?

How far can I go with my ChIP-Seq analysis after Gene Ontology Analysis to find the meaning of my enriched terms and pathways?

Thanks
Benjy

Answer: You can find the up/down regulated genes by performing RNAseq. You're done analysing the ChIPseq data, do the wetlab experiment.
Question: I have a table like below, it is a bed file of genome coordinate, I would like to keep only rows with numbers.

Input:

1	141009669	141009952
9	141016322	141016973
GL000195.1	81719	82468
GL000195.1	142613	142923
GL000220.1	119445	119746
HG115_PATCH 101957832	101958132
HG1308_PATCH 130205069	130205369
HG1308_PATCH 130205406	130205773
HG748_PATCH 77577953	77578264
X	200983	202660
y 205180	205702

output:

1	141009669	141009952
9	141016322	141016973
X	200983	202660
y	205180	205702

Thanks in advance!

Answer: you can do the opposite

awk '{print $1}' file | sort | uniq

to find the rows you want to exclude and then exclude them with grep

grep -v -e 'pattern1' -e 'pattern2'
Question: I download the UCSC genome package named jksrc.zip, how can i install and use it in the Ubuntu linux system?

Answer: http://genomewiki.ucsc.edu/index.php/Source_tree_compilation_on_Debian/Ubuntu