cannin/biostars_qa huggingface dataset:
Question: Hi,

I'm new to Ribo-Seq and trying to understand how the analyze the data. I have difficulties understanding what the offset means. I know it is defined as the distance between the start of the read at the 5'-end and the P-site, which is inside the ribosome.

But how is this offset calculated? why do I have sometimes different offset values for different read lengths, but sometime not and why do i use different offsets for the same read length in different data sets?
Is there a tutorial or a workflow explaining this approach?

thanks

Answer: Hi Assa,

Apologies for making you wait so long for a response. I hope this answer can still be of some use to you.

You are correct with your definition of the Offset for the most part but I believe we can be a small bit more accurate. I think of the offset as the distance from one end of the Ribo-Seq read to the position of the decoding centre of the ribosome on the fragment. Depending on the dataset it may be optimal to calculate offset from the 3' end and the offset value will vary as a result. For elongating ribosomes, such as those isolated when cyclohexamide is used to halt translation we may prefer to calculate offset to the A-site.

There are a number of approaches that have been used to calculate the offset. Historically, set offset lengths were employed for a whole sample depending on the samples treatment. Elongating ribosomes got a set offset of 15 and initiating got a set offset of 12. This was in order to target the P-site for initiating ribosomes and the a-site for elongating

However, this approach is not optimal as the true distance to the decoding centre will differ depending on the length of the fragment. I think Figure S2 in the [ORFik publication][1] is a nice visualisation of this. And a typical calculation of 5' and 3' offset is described in the [Shoelaces paper][2].

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8214792/
[2]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6052522/
Question: I received a script from a colleague and I'm figuring out how samplelist.txt looks like.

cat samplelist.txt |parallel -j 6 --verbose --progress java -jar $(TRIMMOMATIC) PE -threads 4 -phred33 {}_R1.fastq.gz {}_R2.fastq.gz {}_R1.fq.gz {}_U_R1.fq.gz {}_R2.fq.gz {}_U_R2.fq.gz ILLUMINACLIP:trimSeq.fasta:2:7:7 MINLEN:20

These are the fastq files the script is intended to trim:

WT_11-3_R1.fastq.gz
WT_11-3_R2.fastq.gz
WT_12-8_R1.fastq.gz
WT_12-8_R2.fastq.gz
WT_9-30_R1.fastq.gz
WT_9-30_R2.fastq.gz
delta_G_11-3_R1.fastq.gz
delta_G_11-3_R2.fastq.gz
delta_G_12-8_R1.fastq.gz
delta_G_12-8_R2.fastq.gz
delta_G_9-30_R1.fastq.gz
delta_G_9-30_R2.fastq.gz

Answer: In ideal case it should be One sample name per line.

WT_11-3
WT_12-8
WT_9-30
delta_G_11-3
delta_G_12-8
delta_G_9-30
Question: I am wondering could anyone tell me the best work flow (software) and Statistical model should I use for RNA seq data (without replications).

Thanks a lot

Answer: The best tool you'll find is GFold, and you should obviously be very cognizant of the caveats that come with trying to draw any firm conclusions from unreplicated data (i.e., the most this is useful for is getting some target genes for replicated targeted tests).
Question: Dear all,

I am performing an analysis of MEIs in clinical exomes (analogously to https://www.nature.com/articles/s41436-020-0749-x ) and there is something around 0.1% of cases affected by MEIs (as shown in a paper), however, I wander if I miss anything. The questions:

1) so far I detect around 3 different types of MEI using the tool described in the paper ( https://github.com/GeneDx/scramble/blob/master/cluster_analysis/resources/MEI_consensus_seqs.fa ) - are there more? where can I get a catalogue of sequences for different mobile elements?

2) we also detected large retroduplications ( https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3847774/ ) with third-party tools - is there any paper on clinical value of such insertions of whole mRNAs?

3) discussion: which are other potentially genetic disruptions except SNVs/indels/SVs[CNVs]/MEI/retroduplications which may play at least some role in some super rare cases AND are potentially detectable with WES samples? (for cancer we also detect microsatellites - but this is a different story, I don't expect them to play a role in germline)

Thanks a lot!

Answer: Hi,

Thanks for your interest in our paper. We found a diagnostic rate of MEIs in about 0.03% of cases. The 0.1% figure comes from this paper:
https://pubmed.ncbi.nlm.nih.gov/15643617/.

To answer your questions as your enumerated them:

1) We provided fasta sequences for the 3 most active families of MEIs in humans. SCRAMble is sensitive enough to detect other subfamilies based on these sequences (for example, the Alu sequence provided is a consensus AluYa5, but SCRAMble will still be able to detect an AluYb8). If you're interested in using different sequences, I recommend going to RepBase (https://www.girinst.org/repbase/) to find consensus MEI sequences.

2) The DDD cohort has recently published their analysis of MEIs and retroduplications in a clinical cohort (https://www.nature.com/articles/s41467-019-12520-y). They used MELT (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5668948/) to detect both types of insertions. That paper is a good place to start. In our own experience, we've found that rare processed pseudogene insertions affect sensitivity in detecting variants due to dilution of signal from the "real" gene. Similarly, variants that may appear to be of clinical significance might not confirm if they are coming from the processed pseudogene instead of from the original gene.

3) I would add potentially balanced events such as translocations and inversions to this list.

If you have any more questions, you can email me directly by going to the Genetics in Medicine paper and finding the corresponding author info.

Best,
Rebecca
Question: Hello BioStars,

Is there a tool which generates mapping quality score distributions from bam files?

I know I could potentially do this myself, but I am looking for something which would essentially do the work for me, and also something which plugs in nicely to multiqc

I have had a quick search, and it is surprisingly difficult to find a tool which provides this functionality, so I am just wondering if I am missing something obvious here? The obvious candidate would be something like picard, but picard does not seem to offer this functionality from what I can tell

Thanks in advance!

Answer: ```
$ samtools view input.bam | cut -f 5 | sort | uniq -c | sort -n | awk '{printf("MAPQ:%s\t%d\n",$2,$1);}' | gnuplot -e " set terminal dumb ;set nokey; plot '-' using 2:xtic(1) with boxes"

400 ++---+---------+---------+---------+--------+---------+----***********
| + + + + + + * + *
350 ++ * +*
| * *
| * *
300 ++ * +*
| * *
250 ++ * +*
| * *
200 ++ * +*
| * *
| * *
150 ++ * +*
| * *
100 ++ * +*
| * *
| * *
50 ++ *********** +*
| + + + + *********** + * + *
0 **********************************************************************
MAPQ:7 MAPQ:25 MAPQ:0 MAPQ:70 MAPQ:29 MAPQ:37 MAPQ:60

```
Question: I believe Dobzhansky's "Genetics and the Origin of Species" has entered the public domain, or am I wrong? It is from 1937. Likewise for several fundamental books from the Modern Evolutionary Synthesis (Fisher's "The Genetical Theory of Natural Selection", Haldane's "The causes of evolution", etc). However, I can't find eletronic reprints of any of them. Does anyone knows which and where those books are available online?

Answer: (Disclaimer: I'm not a lawyer)

In most countries, copyright is based on the year of death of the author, not the publication date of the work. To quote Wikipedia:

> In most other countries that are signatories to the Berne Convention, copyright term is based on the life of the author, and extends to 50 or 70 years beyond the death of the author. (See List of countries' copyright length.)

Dobzhansky died in 1975 so it will take until 2045 in most countries until his works are public domain. Fisher died in 1962 so it'll take until 2032.

George Orwell, on the other hand, died in 1950, so his works are in the public domain here in Australia where the limit is 50 years (look at ebooks@Adelaide), but not in the public domain in the US, where the limit is 70 years.

Actually, Fisher's work should be in the public domain in Australia (1962+50=2012), The Genetic Theory of Natural Selection is on archive.org, but only the first 1930 edition, not the 1958 second edition (link via Razib Khan). Maybe no-one bothered to put together a 'nice' edition yet? Edit: No wait, I misunderstood - Australia's 50 years expiration only applies if the author died before 1955 (as in Orwell's case), after that it's 70 years like in the US, so Fisher's work is not in the public domain in Australia.
Question: I'm trying to download CDS sequences for a given genome using Biopython. My script looks like this:

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "c...@gmail.com"
genomeAccessions = ['NC_021353.1', 'NC_020913.1']
handle = Entrez.efetch(db="nucleotide", id=genomeAccessions, rettype="gb")
records = SeqIO.parse(handle, "gb")

for i,record in enumerate(records):

print(len(record.features))

for feature in record.features:
if feature.type == "CDS":
print feature.location
print feature.qualifiers["protein_id"]
print feature.location.extract(record).seq

But using this code I get only one feature (for example for genome NC_021353) even though there are many features http://www.ncbi.nlm.nih.gov/nuccore/NC_021353.

I would be grateful for any suggestion what I'm doing wrong.

[1]: http://www.ncbi.nlm.nih.gov/nuccore/NC_021353

Answer: Your rettype must be "gbwithparts", this should work:

handle = Entrez.efetch(db="nucleotide", id=genomeAccessions, rettype="gbwithparts")

Running your program:

>>>
======================== RESTART: D:/Desktop/test.py ========================
3823
[580:1822](+)
['WP_020447944.1']
ATGGACAACATCAGCGGCAACATCTTTG...
...
[1818:2619](-)
['WP_020447945.1']
GTGAATAAAGAGAACAAATCTACTTTAA...
...

Your code will break later, possibly there are entries without a protein_id:

Traceback (most recent call last):
File "D:/Desktop/test.py", line 20, in
print feature.qualifiers["protein_id"]
KeyError: 'protein_id'
Question: Hello everyone,

I have two dataframes: df1 and df2. I want to combine them using `rbind` function. These two dataframe have some overlapping rownames. When combining it adds a value to the overlapping rownames.

For example:

df1 has one row as rs119151

df2 also has one row as rs119151

After combining its coming as

rs119151
rs1191511

I want it to show as:

rs119151
rs119151.1

and so on if we have more same rows for this value.

I am trying to use make.names command but its not working for me?

```r
new <- make.names(rbind(df1,df2),unique=TRUE)
```

Is there any other way to do it?

Answer: `make.names` works on a vector, not a dataframe.

We can *rbind*, then reassign fixed *rownames*, try this example:

d1 <- data.frame(x = 1:2, row.names = c("rs1", "rs2"))
d2 <- data.frame(x = 4:5, row.names = c("rs1", "rs3"))

res <- rbind(d1, d2)
# x
# rs1 1
# rs2 2
# rs11 4
# rs3 5
rownames(res) <- make.unique(c(rownames(d1), rownames(d2)))
res
# x
# rs1 1
# rs2 2
# rs1.1 4
# rs3 5
Question: Hi,
I want to make a data table that I have, to something that is more comprehensive. In abstract this is what I intend to do:

> df
cond_1 cond_2 cond_3 cond_4
1 a a b e
2 b b c f
3 c d g

> df_t
cond_1 cond_2 cond_3 cond_4
a + + - -
b + + + -
c - + + -
d - - + -
e - - - +
f - - - +
g - - - +

Is there any package or function that can do the job or facilitate coding? Any suggestion is a big help!

Answer: rownames = unique(unlist(df))
rownames = rownames[!is.na(rownames)]
df2 = sapply(df, function(x) ifelse(rownames %in% x, "+", "-"))
df3 = as.data.frame(df2, stringsAsFactors = F, row.names = rownames)

> df3
cond_1 cond_2 cond_3 cond_4
a + + - -
b + + + -
c - + + -
d - - + -
e - - - +
f - - - +
g - - - +
Question: Hi everyone

I need to understand Minor allele frequency

A SNP with a allele (G) frequency of 0.40 implies that 40% of a population has G allele versus the most common allele (major allele), which is found in 60% of the population.

If a SNP allele frequency is 0.05 how will be calculate its percentage in the population

Please explain with calculation

Regards

Answer: If a SNP allele frequency is 0.05 in a population of 100 people:

Total number of alleles for each SNP: 100 * 2 = 200 (each individual has two alleles)

Total number of alleles for your SNP with a MAF of 0.05: 200 * 0.05 = 10

The minor allele is present ten times in your population. However, you don't know how many heterozygote and homozygote.
Question: Hi Everyone,

I want to get geneID(gene_id) and its associated gene symbol(gene_name) from a gtf file into a file that can be read into R.

My gtf file looks like this:

```
chr1 HAVANA gene 3073253 3074322 . + . gene_id "ENSMUSG00000102693.1"; transcript_id "ENSMUSG00000102693.1"; gene_type "TEC"; gene_status "KNOWN"; gene_name "RP23-271O17.1"; transcript_type "TEC"; transcript_status "KNOWN"; transcript_name "RP23-271O17.1"; level 2; havana_gene "OTTMUSG00000049935.1";
chr1 HAVANA transcript 3073253 3074322 . + . gene_id "ENSMUSG00000102693.1"; transcript_id "ENSMUST00000193812.1"; gene_type "TEC"; gene_status "KNOWN"; gene_name "RP23-271O17.1"; transcript_type "TEC"; transcript_status "KNOWN"; transcript_name "RP23-271O17.1-001"; level 2; tag "basic"; havana_gene "OTTMUSG00000049935.1"; havana_transcript "OTTMUST00000127109.1";
```

Can anyone tell me how to write a little script (preferably python) to do it?

Any idea would be appreciated!

Answer: One quick and dirty way to extract attributes:

```py
#!/usr/bin/env python

import sys

for line in sys.stdin:
attr = dict(item.strip().split(' ') for item in line.split('\t')[8].strip('\n').split('') if item)
sys.stdout.write("%s\n" % (attr['gene_id'].strip('\"') + '\t' + attr['gene_name'].strip('\"')))
```

To use:

$ ./extract_gtf_attributes_of_interest.py < foo.gtf > answer.txt

The result with your test input:

```
$ more answer.txt
​ENSMUSG00000102693.1 RP23-271O17.1
ENSMUSG00000102693.1 RP23-271O17.1
```

As mentioned, this is quick and dirty and does no error checking or input validation. Consider adding a `try..except` block to handle the case where an attribute field is missing a value (or the attribute field is missing altogether).
Question: I realized that I don't know what happens when illumina sequencing chemistry reaches the end of a fragment. Does the reaction stop for that fragment, or are bases added in some way? The reason I ask is because I have fragments ranging from 120-160bp in length and yet on 150 cycles I will always get 150bp long reads. And some of my reads will end in terminal repeats like this:

CGTCTTCTGCTTGAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

Answer: "CGTCTTCTGCTTG" matches the end of an Illumina adapter sequence. After you run off the end of an adapter, it's common to get poly-A, and sometimes eventually poly-G. So, no, your adapters have NOT been trimmed - you are sequencing into them, and then off the end into no signal.

You can (and should) trim your adapters, as described in [this post][1].

Also, to clarify, the Illumina sequencing machine has no idea how long your fragments are. You tell it to run for 150 cycles, so it gives you 150 cycles of data for every read, regardless of the actual length of the molecule.

[1]: http://seqanswers.com/forums/showthread.php?t=42776
Question: ~~I am trying to update my bioconductor to the latest version (3.0+). I currently have R 3.2.1 and Bioconductor 2.28. When I run BiocUpdate, it just reinstalls 2.28. I tried uninstalling and reinstalling Bioconductor, but still only got 2.28.~~

~~The reasoning is that a package I am running (phyloseq) is having an error that the developer doesn't recognize. He suggested updating to the latest version of Bioconductor & R to fix it. I have updated my R version, but am unable to update my Bioconductor.~~

~~Any suggestions to install Bioconductor 3.0 or 3.1?~~

**Edit: There was a misunderstanding in how to tell what version of Bioconductor I have.**

**The best way to check the version of Bioconductor is to compare your version of BiocInstaller (a package within R once Bioconductor is installed) with the Bioconductor website.**

Answer: The best way to check the version of Bioconductor is to compare your version of BiocInstaller (a package within R once Bioconductor is installed) with the Bioconductor website.
Question: I have a number of VCF files, where each VCF file possesses variant data for a single patient (this is the way Illumina provides their data). Is it possible to combine all of the data for the patients into one VCF file? If so, how? Can I use plink/seq to do this?!

Any suggestions and leads would be extremely helpful.

Answer: Related duplicate post

- https://www.biostars.org/p/49730/

Use [vcf-merge](https://vcftools.github.io/perl_module.html#vcf-merge)
Question: I have met an weird error.

Here lists my two bed files

`a.bed`

```
chr1 645339590 645344590
chr1 645485603 645490603
chr1 646267584 646272584
chr1 651156159 651161159
chr1 652224601 652234601
chr1 654584466 654589466
chr1 656200028 656210028
chr1 658335720 658345720
chr1 672881471 672886471
chr1 673481036 673486036
```

`b.bed`

```
chr1 567457 567702
chr1 569801 570047
chr1 713857 714086
chr1 714475 714736
chr1 762251 763163
chr1 762416 762594
chr1 762691 762997
chr1 793349 794537
chr1 800991 801931
chr1 805262 805484
```

When I was running `intersectBed -a a.bed -b b.bed`, there is an error `intersectBed: line 2: 21275 Segmentation fault (core dumped) ${0%/*}/bedtools intersect "$@"`.

When I was running `intersectBed -a b.bed -b a.bed`, there is another error `ERROR: Received illegal bin number 39388 from getBin call. ERROR: Unable to add record to tree.`

I can not figure out what causes this. Thanks for the help.

Chen Tong

Answer: @aaronquinlan will confirm but I suspect that the value for chr1 is just to high for the internal 'binning' system of bedtools. See http://genomewiki.ucsc.edu/index.php/Bin_indexing_system

position max = 536,870,912

but you used

position max = 673,486,036

my idea / hypothesis.
Question: I've got some data with p values expressed as -log(p). How do I invert this to the p values?

I've read here about undoing -log10(p):

https://www.biostars.org/p/358046/

```r
p <- 0.01
logp <- -log10(p)
10^-logp
```

Here are some example -log(p) in the table:

```
1.542108761
1.760644949
1.567865477
1.779197974
```

Is it:

exp(1)^(-logp)

Thanks!

Answer: In base 10, it should be:

```r
1/(10^num1)
```

From example:

```r
p = .05
num1 = -log10(p)
1/(10^num1)
[1] 0.05
```
Question: Hi,

I have a multi-fasta file with a header in the following format:

>sp|Q9Y5Q8|TF3C5_HUMAN General transcription factor 3C polypeptide 5 OS=Homo sapiens GN=GTF3C5 PE=1 SV=2

I would like to format to extract the Uniprot ID or the Protein Name (ACC) to get the following:

>Q9Y5Q8

or

>TF3C5_HUMAN

I think sed can do it but I don't know the exact combination of regexp

Thanks

Answer: awk '{if ($0 ~ /^>/) {split($0,a,"|"); print ">"a[2]} else { print;}}' your_file > new_file

If you want the `TF*` names then

awk '{if ($0 ~ /^>/) {split($0,a,"|"); split(a[3],b," "); print ">"b[1]} else { print;}}' your_file > new_file
Question: I have HLA typed whole genome sequencing data using [HLA-VBSeq][1]. Reads were first aligned to hg19 with BWA-MEM, then reads aligning to all HLA loci and unmapped reads were extracted from the bam file using samtools. The extracted reads were then combined and re-aligned to the collection of all the genomic HLA allele sequences in the IMGT/HLA database using BWA-MEM, allowing for multiple alignments to the reference sequences.

Now, I would like to subset the reads of a specific HLA-E gene allele, however, using `samtools view input.bam HLA:HLA00936`, I see that all of the sequences have a MAPQ score of 0. I know that this is specific to BWA and indicates that the reads are mapping to multiple locations. Does this mean that the alignment results for this allele are insignificant and should be disregarded? The researcher interested in this allele would like to know the HLA-E sequences for gene editing purposes, but I feel as though the coverage is far too low. Any input or advice is much appreciated.

[1]: http://nagasakilab.csml.org/hla/

Answer: If you have multiple inferred HLA-E alleles, then it is expected to have most reads mapping to multiple sequences, because HLA-E is one of the least polymorphic HLA genes. Many sequence segments will be identical between the alleles. You can get the [nucleotide sequences of the HLA-E alleles][1] as a FASTA file and see how similar they are with the multiple sequence alignment tool [MUSCLE][2] (choose HTML output format). You must have some reads mapping uniquely to the allele you are investigating because you inferred that the allele exists from the data.

[1]: ftp://ftp.ebi.ac.uk/pub/databases/ipd/imgt/hla/fasta/E_nuc.fasta
[2]: http://www.ebi.ac.uk/Tools/msa/muscle/
Question: I have a long list of Ensembl IDs and am completely unfamiliar with the Ensembl browser. I want to convert these guys to Refseq IDs using the UCSC table browser, mySQL, or some other tool. So given a list like this:

ENSG00000197021
ENSG00000204379

How do I turn them into this?

CXorf40B
XAGE1A

I gather from what I've tried that there are not 1-to-1 mappings for Ensembl gene IDs to RefSeq names, but I'd like to be wrong. At a minimum, I'd like to get the hg18 or hg19 coordinates for a list of these genes, but Biomart only comes back with about 3/4 of them. Any help?

Answer: using the public mysql server of UCSC:

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19
mysql> select distinct G.gene,N.value from ensGtp as G, ensemblToGeneName as N where
G.transcript=N.name and
G.gene in ("ENSG00000197021", "ENSG00000204379") ;
+-----------------+----------+
| gene | value |
+-----------------+----------+
| ENSG00000197021 | CXorf40B |
| ENSG00000204379 | XAGE1A |
+-----------------+----------+
2 rows in set (0.21 sec)
Question: I am wondering what the best practices are analyzing host-parasite mixed transcriptomes. If I have mRNA-Seq data from infected host tissues that contain both host and parasite genomes, which would be the best way to tease them apart? Mapping to a combined reference (host + parasite), or mapping to the host genome first, collecting unmapped reads and mapping to parasite genome? What would be the implications of these options? When I get the real data, I could run tests with both approaches, but maybe someone already have some thoughts on this here?

Answer: In my case, I mapped the RNA-seq reads to both host and parasite genomes simultaneously, after concatenating their genomes. Then, I removed unmapped reads and reads with primary alignment to the host chromosomes (I was interested in the parasite). From my point of view, if you remove reads that map to the host genome without verifying if they map to the parasite genome, you could end up removing reads derived from parasite's genes that have highly conserved homologs in the host genome, even though such cases are rare.
Question: Hi,

I' working on the protein (cytochrome oxidase), that is a complex of 3 to 13 chains. There is a very good model 3D model. The first subunit is quite big and has few transmembrane domains. I would like to choose one amino acid and check which other amino acids are close to my chosen one (there is an interaction between them). Which software could show that?

Thanks,

Answer: Pymol can do it, as can Chimera or VMD or other visualization program. You can also script it if visualization is not necessary. [In Pymol, you load your structure and use the "around" or "within" selectors][1].

For example, to select residues within 5Å of histidine 34 you would write the following:

select my_selection, byres( resi 34 around 5 )

[1]: http://www.pymolwiki.org/index.php/Selection_Algebra