cannin/biostars_qa huggingface dataset:
Question: Hello community!

I just started with bioinformatics, so I hope you can help me with understanding some basic things. I have RNA-Seq 100bp pair-ended reads from the *Lepidoptera *insect, for which there is no genome available. In total there were 2 treatments and one control, so for the assembly I merged all of them into two files (R1 and R2). So the goal is to make an assembly, and then continue with mapping, counts and statistic.

For the assembly I used next tools: Trinity, Velvet/Oases, Bridger and CLC assembler. From the last one I just got the ready-to-go assembly, to which I compared other assemblers.

1. Trinity was executed with default parameters (k=25),
2. Velvet: velveth/velvetg k=21,51,2; For each assembly run oases and compared n50 and total number of contigs. Chose the range of contigs from 25 till 39, with kmer=35, velveth/velvetg on the oases results and then merged them with oases.
3. Bridger. Default parameters with kmers 25 and 27.

After I compared assemblers with Quast. As a result got the next table

```
Oases Trinity Bridger 25 Bridger 27 CLC
#contigs 10706 15579 14881 14632 19985
#contigs (>= 0 bp) 36576 40429 37011 36162 35589
Largest contig 27732 18062 27729 27887 40347
Total length 12451630 19701362 20030913 19559912 36090599
N50 1324 1528 1653 1651 2584
# N's per 100 kbp 0 0 0.05 0.21 212.61
# N's 0 0 11 41 76733
```

So, I could not reach the N50 of CLC, by using default parameters. But at least I can see how do they reach it. The number of mismatches in CLC is super high, so I was wondering - is it really worth to increase the total length and n50 by putting double amount of N's into the contigs?

And the next question - for some reason I though it would be a good idea to increase the input sequence length by merging R1 and R1, R2 and R2 reads. So basically doubling the amount of data assembler should deal with. So far I tested only Bridger (k=25), but the results are surprising, N50 increased to 1720, and the total length and the largest contigs became almost the same as in CLC results. What I think is happened that during the building of the graph low quality k-mers had a better chance to survive, thus increasing the contig length.

Thanks in advance for any kind of responds!

Answer: Mapping back the reads is definitely one of the good way to evaluate the assembly. Apart from it, you can try:

1. Comparing the assemblies by performing [BLAST][1] them against nearest species full-length cdna's / proteins sequences obtained from RefSeq and check how many your transcriptome picked up.
2. Finding full-length sequences in your assembly (using something like [full-lengther_next][2])
3. Finding Open Reading Frames in each of your assemblies.

[1]: http://blast.ncbi.nlm.nih.gov/Blast.cgi
[2]: http://www.rubydoc.info/gems/full_lengther_next/frames
Question: Hello,

I have two questions regarding merging the `BAM` files. I have several libraries 4 lanes each. I demultiplexed, mapped against the reference genome, gave the read group name. I am using `snakemake` and `samtools merge`.

These are my input files:

Bob_1_CL1.bam
Bob_2_CL1.bam
Bob_3_CL1.bam
Bob_4_CL1.bam
John_1_CL40.bam
John_2_CL40.bam
John_3_CL40.bam
John_4_CL40.bam

This is the code I am trying:

library=["Bob", "John"]
(LIBRARIES, SAMPLES)=glob_wildcards("path/to/{library}_1_{sample}.bam")

rule all:
input:
expand("{library}_{sample}.merged.bam", library=LIBRARIES, sample=SAMPLES)

rule samtools merge:
input:
lane1="path/to/{library}_1_{sample}.bam",
lane2="path/to/{library}_2_{sample}.bam",
lane3="path/to/{library}_3_{sample}.bam",
lane4="path/to/{library}_4_{sample}.bam"
output:
outf = "{library}_{sample}.merged.bam"
threads: 4
shell:
"""/Tools/samtools-1.10/samtools merge -@ {threads} {input.lane1} {input.lane2} {input.lane3} {input.lane4} -o {output.outf}"""

I am trying to merge 4 lanes in one single `BAM` file in order to get `Bob_CL1.merged.bam` and `John_CL40.merged.bam`, the problem is that CL number (NOTE: it is not always CL) is different for different libraries, i.e. CL1 is present in Bob, but not in John, and CL40 is for John but not in Bob. so it is giving me an error:

Missing input files for rule merge:
John_1_CL1.bam
John_2_CL1.bam
John_3_CL1.bam
John_4_CL1.bam
So my questions are:

1. How to solve this?
2. Do I have to merge the lanes in one single `BAM` file prior to marking duplicates?

Answer: 1. If you print the output of `expand("{library}_{sample}.merged.bam", library=LIBRARIES, sample=SAMPLES)` you'll see the problem. Replace `expand()` with something like `['{}_{}.merged.bam'.format(x, y) for x, y in zip(LIBRARIES, SAMPLES)]`
2. Yes, PCR duplicates are library-level events, so you need to process the entire library at once.
Question: Hi,

I have an alignment between 4 proteins and I want to add each protein's sequence homolgs to it. Is there a tool where I can label my original alignment as a "core alignment" and not change it, and have the homolgs align according to this alignment?

Thanks!

Answer: MAFFT lets you do this with the `--add_sequence` option: http://mafft.cbrc.jp/alignment/software/addsequences.html

Quite good alignment software, and I've used this option many times to add new sequences to existing alignments.
Question: With this input parameter:

fci__rf: [FailsVendorQualityCheck, BadMate, UnmappedRead, BadCigar]

I would like this:

read_filters:
type: string[]
inputBinding:
prefix: --read_filter

To result in this:

java
-jar GenomeAnalysisTK.jar
-T
FindCoveredIntervals
--input_file
sample.bam
--reference_sequence ref.fasta
--read_filter FailsVendorQualityCheck
--read_filter BadMate
--read_filter UnmappedRead
--read_filter BadCigar

Do I need to use `shellQuote`? Or perhaps I need to use a string input instead of an array?

Answer: Please check out this page http://www.commonwl.org/user_guide/09-array-inputs/

In your case this definition will do the job:

inputs:
read_filters:
type:
type: array
items: string
inputBinding:
prefix: --read_filter
inputBinding:
position: 1
Question: I am a biologist do not have very good exposure to bioinformatics.

Recently I came across some of my colleagues are using Galaxy Australia to analyze RNASeq data.

I had analyzed RNASeq data using R Programming Bioconductor packages (EdgeR, and DESEQ) along with FastQC, Trim_Galore, Kallisto. It was challenging to get access to servers and learn Linux.

May I know what is Galaxy Australia and how it has advantages over other pipelines?

I am interested in doing an RNASeq for a few crops which are not well-annotated genome.

Answer: I had analyzed RNASeq data using R Programming Bioconductor packages (EdgeR, and DESEQ) along with FastQC, Trim_Galore, Kallisto. It was challenging to get access to servers and learn Linux.

Most of these tools are available on Galaxy (probably on Australia Galaxy) and you can drag and drop the tools to construct a workflow. Once workflow is built, customize the parameters on UI and execute the workflow

May I know what is Galaxy Australia

Galaxy is a web based, click and run data analysis platform to execute tools and workflows. Galaxy Australia is a server hosting Galaxy platform, maintained by [this team](https://usegalaxy-au.github.io/people) based on Australia.

and how it has advantages over other pipelines?

Galaxy Australia or Galaxy it self is not a pipeline. It's a data analysis framework. Compared to other frameworks (esp NGS data analysis), learning curve is minimal with Galaxy server for creating, editing and executing workflows. It's advantages are drag and drop building of workflows, visual editing of workflows, sharing with users, export the workflow and history. In addition, user can search, import and execute workflows shared by community members within Galaxy.
Question: Hello everyone!

I'm a newbie in analyzes with RNA-Seq.

I have paired-end reads data (Sanger / Illumina 1.9) and would like to cut the reads using quality score.

We have visualized data with FastQC and we tested the softwares: sickle, Spade, fastq quality filter (fastx) and AllPaths.

The best result was obtained with the fastq quality filter, but we lost a lot of reads in the process (it does not cut off part of the reads, only excludes the entire read).

I wonder if there is already a program that recognizes the first occurrence of a symbol (eg '#') and cuts the sequence and the quality from that point.

Answer: When a major part of the read has low quality bases then trimming reduces the length of the read and now the read can't be aligned with higher confidence against the reference genome. So almost all the trimming software will discard reads whose length has been reduced to less than some number (say n=30). If most of your reads suffer from this problem than all the trimming tools will behave in the same way i.e. discarding majority of the reads. Although this is not a solution but I would suggest you trying [Trimmomatic][1] and see if it helps.

[1]: http://www.usadellab.org/cms/?page=trimmomatic
Question: Hey guys, I have a table with circa of 1600 rows, like this:

> Element-ID,Protein-ID,Protein-Product,Virus-ID,Super-Kingdom,Order,Family,Genus,Species,Sense,Start,End,EVE-length
> ctg_1003:180657-180956,YP_184770.1,hypothetical protein,39640,Viruses,Unclass,Polydnaviridae,Bracovirus,Unclass,neg,180657,180956,299
> ctg_1003:283549-284079,YP_184879.1,hypothetical protein,39640,Viruses,Unclass,Polydnaviridae,Bracovirus,Unclass,pos,283549,284079,530
> ctg_1007:58711-59043,YP_184770.1,hypothetical protein,39640,Viruses,Unclass,Polydnaviridae,Bracovirus,Unclass,neg,58711,59043,332
> ctg_100:908810-909199,YP_184882.1,hypothetical protein,39640,Viruses,Unclass,Polydnaviridae,Bracovirus,Unclass,neg,908810,909199,389
> ctg_1011:242875-243240,YP_001426207.1,hypothetical protein,399781,Viruses,Algavirales,Phycodnaviridae,Chlorovirus,Paramecium,bursaria Chlorella virus A1,pos,242875,243240,365

Basically it's a csv file with taxonomy information from endogenous viruses. So I read this csv file as a pandas dataframe (I'm working with python 3.7).

df = pd.read_csv('eve_tax.csv', sep=',',header = 0)

An create a swarmplot with seaborn:

beeswarn_plot =sns.swarmplot(x='EVE-length',y='Super-Kingdom', hue='Family', data=df)

So, my code, at the moment, is:

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.style as style

sns.set(style="ticks")
style.use('seaborn-poster')

df = pd.read_csv('eve_tax.csv', sep=',',header = 0)
beeswarn_plot =sns.swarmplot(x='EVE-length',y='Super-Kingdom', hue='Family', data=df)
sns.despine(fig=None, top=True, right=True, left=True, bottom=False, offset=None, trim=False)
beeswarn_plot.set_ylabel('')
beeswarn_plot.set_yticks([])
beeswarn_plot.set_xlabel('EVE length (pb)')

beeswarn_plot.legend(loc='upper center', bbox_to_anchor=(0.7, 1.05), ncol=4, fancybox=True, prop={'size': 11},title='EVEs Family')
plt.savefig('eves_tax_beeswarn_plot.pdf',dpi=300)

Which plots a 'beeswarm' plot like this:

![beeswarmplo][1]

[1]: https://i.ibb.co/sV0d9Vf/Screenshot-from-2020-06-30-08-42-29.png

So, what I want?

For some families I have a low number of elements, such as: Orthomyxoviridae n = 5 , Poxviridae n= 5. I want to put all families with 5 or less elements in a category 'Others', and plot with seaborn all families with 5 or more elements and the category 'Others'.

Can anyone help? Thanks !

Answer: I resolve this with a quick fix:

family_count = df['Family'].value_counts()
family_low = family_count[family_count <= 5]
family_low_list = family_low.index
for i in family_low_list:
df['Family'] = df['Family'].replace({i:'Others'})

So I replace all families that occurs less than 5 times in dataframe for 'Others', but I know that can have a way to make this without quick fixes...
Question: Hi,

This appears to be a simple problem that I am unable to solve. I have some data that looks like this:

CHROM POS REF ALT TYPE AF
chr1 1 A T MISSENSE 0.23
chr2 1 A T,G MISSENSE 0.17, 0.09

The above is dummy meaningless data, but it is representative of the problem at hand.

I'd like to `separate_rows` such that the `ALT` and `AF` are separated in a couples manner. Running `separate_rows` on the 2 columns would give me 4 rows, not 2. I'd like my output to be:

CHROM POS REF ALT TYPE AF
chr1 1 A T MISSENSE 0.23
chr2 1 A T MISSENSE 0.17
chr2 1 A G MISSENSE 0.09

Is there any way I can conserve this combination while separating the values out? I am really far out from the VCF to go back and split multi-allelics.

Answer: If I run this:

separate_rows(data, ALT, AF, convert = TRUE)

where `data` is your 1st data frame I obtain the second data frame.

How did you run the function `separate_rows()`?
Question: Colleagues,

I have some BAM files with excessive coverage in certain regions that I want to cap to a given coverage. In regions with lower coverage, I do not want the coverage to change. Here's an illustration of what I want:

I want read pairs to be intact in the output file.

Does anyone have a solution to this? I want to output a BAM file for other tools do use. Google has not been my friend.

(note: GATK PrintReads using `-dcov` looked promising until a realised that "For read-based traversals (ReadWalkers like BaseRecalibrator), it controls the maximum number of reads sharing the same alignment start position" from https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_engine_CommandLineGATK.php#--downsample_to_coverage. Close but no CIGAR (pun intended)).

Answer: I wrote two tools:

- The first sorts a BAM on REF/read-name https://github.com/lindenb/jvarkit/wiki/SortSamRefName
- The second reads the BAM above. For each set of reads having the same name, it checks if after adding those reads, the depth would be lower than the expected depth. If true, those reads are added to the BAM. **Edit**: doc is https://github.com/lindenb/jvarkit/wiki/Biostar154220

```
$ java -jar dist/sortsamrefname.jar -T bam /input.bam |\
java -jar dist/biostar154220.jar -n 20 -T bam |\
samtools sort - output
```

```
$ samtools mpileup output.bam | cut -f 4 | sort | uniq -c
12692 0
596893 1
94956 10
56715 11
76947 12
57912 13
66585 14
51961 15
63184 16
47360 17
65189 18
65014 19
364524 2
169064 20
72078 3
118288 4
54802 5
82555 6
53175 7
78474 8
54052 9
```

~~(what are those 21 and 22 ? I suppose it's because I don't consider supplementary, secondary, duplicate reads )~~ *(bug fixed)*. Anyway, can you tell me if that's what you need?

Pierre
Question: Hi all,

I am very knew to all things Chip-seq. We have performed multiple experiments and now I have to analyze multiple files and was told to pretty much merge the replicates and use the mean of the distance for each gene.

Up until here I agreed with, but while I understand merging the replicates I was also told to merge broad and narrow peak files.

For all I've been reading this sounds like a terrible idea, but I'm the new person. If I'm correct (there is always the chance that I'm completely wrong) what arguments would be best what could I use as reference to support this position ?

Thank you

Answer: First, are these A) technical or B) biological replicates? That is, the same biological sample run several times with the same antibody (same lot also if polyclonal) protocol, or different biological samples run the same way with the same protocol?

If it is A it may be reasonable to merge them for some analyses, such as just annotating peaks. I would merge the bam alignment files and then do the calls versus merging the calls.

However, first you have analyze your replicates to check they they all perform the same. We did a lot of performance comparisons here:
https://epigeneticsandchromatin.biomedcentral.com/articles/10.1186/s13072-016-0100-6

You can steal those ideas, especially using the ENCODE segmentation tracks if it's human and they have tracks for something like your cell type. But just counting the reads in bins and then doing a correlation is pretty informative.

But even in our data, and we used a robot and do it a lot, one of our technical replicates behaved strangely. See supplemental figure S6.

If it is B, biological replicates, you almost certainly don't want to merge them. You will lose your information about biological variance is present. If you are looking at something like differential peaks between conditions DESeq and really all reputable programs will want some sort of replicates, almost always biological. In general, if you want to compute a p value on anything you need separate replicates (not merged).

If you are just annotating peaks you don't need a p value.
Question: Hello. I've been recently trying to programatically convert a bunch of UniProt IDs to gene names.
I found the UniProt API, which should do the job, something in the lines of:

import urllib,urllib2

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

params = {
'from':'ACC',
'to':'P_REFSEQ_AC',
'format':'tab',
'query':'P13368 P20806 Q9UM73 P97793 Q17192'
}

data = urllib.urlencode(params)
request = urllib2.Request(url, data)
contact = "" # Please set your email address here to help us debug in case of problems.
request.add_header('User-Agent', 'Python %s' % contact)
response = urllib2.urlopen(request)
page = response.read(200000)

The problem is, this returns whole website. Is it possible to only obtain a e.g. JSON where a list of mappings and corresponding information would be present (e.g. species too).

Thank you.

Answer: You can use requests like

http://www.uniprot.org/uniprot/?query=accession:P13368&format=tab&columns=genes

to access only the gene names. However, the delimiting is a bit odd to parse in this case. i.e. some entries are linked to more than one gene. And genes often have more than one name and its hard to figure whats what in this output.

You can either parse this out of the different file formats or use our [sparql][1] endpoint to just ask the preferred gene names
directly.

BASE
PREFIX skos:
PREFIX up:
PREFIX taxon:
SELECT ?protein ?preferredGeneName
WHERE
{
VALUES ?protein { }
?protein a up:Protein ;
up:encodedBy/skos:prefLabel ?preferredGeneName .
}

You can use the download links for this to get the information back as json/xml or csv as you wish and by editing the UniProt accessions in the query you can retrieve all entries you want.

[1]: http://sparql.uniprot.org/
Question: I'm trying to access identical protein information for a list of protein accession numbers. My goal it to take protein ID's, get identical protein list, use the nucleotide coordinates from that list to specifically download the upstream and downstream region of a given gene.

I was hoping to use the biopython Entrez module.

I've used the below code

epost_1 = Entrez.read(Entrez.epost("protein", id=",".join(Prot_ID_list)))
webenv = epost_1["WebEnv"]
query_key = epost_1["QueryKey"]
prot_link = Entrez.read(Entrez.elink(dbfrom="protein", db="protein", LinkName="protein_protein_identical", webenv=epost_1["WebEnv"], query_key=epost_1["QueryKey"]))

However, this doesn't work. It just returns me a list of protein IDs that I originally posted.

I've also tried using eutils on ubuntu

epost -db protein -id Prot_ID_List | elink -related -name protein_protein_identical | efetch -format text

but it just gives an error "QueryKey value not found in fetch input". I tried variations on this to no avail. I am probably doing something very wrong with eutils, but I am not sure what.

Any help is greatly appreciated.

Answer: It appears I was approaching the question wrong. You don't need to 'link' to the identical protein list, despite the 'protein_protein_identical' link name, and description [here][1] on the NCBI website.

Instead, it is information intrinsic to each accession number.

epost_1 = Entrez.read(Entrez.epost("protein", id=",".join(prot_list)))
webenv = epost_1["WebEnv"]
query_key = epost_1["QueryKey"]
iden_prots = Entrez.efetch(db="protein", rettype='ipg', retmode='text', webenv=epost_1["WebEnv"], query_key=epost_1["QueryKey"])

I used the above to get the identical protein list. rettype=ipg downloads the identical protein list in tab-delimited format.
[1]: https://eutils.ncbi.nlm.nih.gov/corehtml/query/static/entrezlinks.html

Hopefully other people come across this if they get stuck like I did.
Question: After checking exons length from the gencode.v15.annotation I have noticed that there are exons only 1bp or 2bp in length.

curl -s "ftp://ftp.sanger.ac.uk/pub/gencode/release_15/gencode.v15.annotation.gtf.gz" |
gunzip -c |
awk '($3=="exon" && $5-$4+1 < 3) {print}'

Thats strange as some exons (protein coding or lncRNA are only 1bp or 2bp long). Is it bioinformatics or (probably not) biology? Has anyone ever noticed something like that with different annotation?

Answer: I have contacted and asked Gencode staff about this issue (in February). They have answered and hoped that problem will be fixed until Gencode.v16.

Apparently there was a bug in one of their scripts.

"... there should be no exons in Gencode <3bp. Alignments of <3bp can not be trusted, even when spanning known splice junctions, or confirming known UTRs/retained introns".

Current Gencode annotation (v18) still have this problem (don't know why they haven't fixed it yet).

I would suggest filtering those exons out.
Question: I am using autolplot function from the ggfortify library in R. Autoplot serves cluster analysis too. I wanted to know what is the algorithm used by the autoplot for finding the 1st two principal components?

Answer: The PCA function that it uses is `prcomp()`, which is the same as what my own package (PCAtools) and DESeq2 use.

Yes, it is performing partitioning around medoids (PAM) and identifying X number of clusters (user pre-selects desired number as second parameter to `pam()`). `autoplot()` then performs PCA on the dataset and shades the points based on the PAM cluster assignments. Here is the proof:

g1 <- autoplot(prcomp(iris[-5]), frame = TRUE, frame.type = 'norm')
g2 <- autoplot(pam(iris[-5], 2), frame = TRUE, frame.type = 'norm')
require(grid)
require(gridExtra)
grid.arrange(g1,g2, ncol = 2)

They are the same points, but higlighted differently.

As is typical with many CRAN (and other) packages, the documentation is poor and the program functionality does not make it readily obvious what the function is doing.
Question: Hello again,

I was curious if anyone had a nice way of finding the gene synonyms for a list of genes? My genes are currently listed as gene symbols, however, any tool or method to fetch synonyms from any type of ID would be a step in the right direction. Conversions in the reverse direction would also be helpful (from any known synonym to gene symbol).

As an example: If I came across IL8 (a synonym), I would convert it to CXCL8 (the gene symbol), or vice-versa. I know in ENSEMBL (and many other tools) you can search for genes using many synonymous names, so there must be a table somewhere to easily covert between the two, but I have not been able to locate it.

Currently, I am using a massive table from UniProt and using 'awk' to search for the genes I need. This is both cumbersome and inefficient, which is why I wanted to reach out to you all!

Thanks!

Answer: using mysql / ucsc

~$ mysql --user=genome --host=genome-mysql.soe.ucsc.edu -A -P 3306 -D hg38 -e 'select distinct K2.alias from kgAlias as K1, kgAlias as K2 where K1.alias="IL8" and K1.kgID=K2.kgID;'
+-------------------+
| alias |
+-------------------+
| B2R4L8 |
| CXCL8 |
| ENST00000307407.7 |
| IL8 |
| IL8_HUMAN |
| NM_000584 |
| P10145 |
| Q6FGF6 |
| Q6LAE6 |
| Q96RG6 |
| Q9C077 |
| Q9UCE1 |
| Q9UCR8 |
| Q9UCR9 |
| Q9UCS0 |
| uc003hhe.1 |
| uc003hhe.2 |
| uc003hhe.3 |
| C9J4T6 |
| C9J4T6_HUMAN |
| ENST00000401931.1 |
| LP896310 |
| uc062xgp.1 |
+-------------------+