cannin/biostars_qa huggingface dataset:
Question: Hi all,

I recently downloaded HMMER to use **hmmscan** locally in Linux with a Pfam database. It works great, however the ouput files are quite difficult to read quickly in my opinion...

I tried using output options such as: --tblout, --domtblout, --pfamtblout, etc. but the ouput files are still voluminous.

**I would like to keep only the top hits in my output files.**

I've seen that it was possible with hmmsearch so I was wondering if there was something similar with hmmscan...
Ideally, I would want an output as I could find online: ![cf. here][1]

If you have any suggestions, I'll gladly took them. Thank you in advance for your very appreciated help!

[1]: https://iili.io/d0LcrX.png

Answer: For anyone interested, I figured it out using this amazing resource: **http://slhogle.github.io/2015/remove-duplicate-lines/** and the option **--tblout** of hmmscan.

I did:

hmmscan --tblout output_file.pfam Pfam-A.hmm seq_file.fasta

and then:

awk '!x[$3]++' ouput_file.pfam > MYBESTHITS.pfam

MYBESTHITS.pfam file is basically what I got online with the top hit for each protein sequences.
Question: Dear group,
Given RNA-Seq data for either 1 or more than 1 different samples/conditions, I am interested in finding

A. an exon-intron or intron-exon junction with reads spanning these junctions.
B. Another similar question is identify intron-retention, where reads map exon-intron-exon.

A differs from B, where A coverage ends halfway in intron. There are algorithms such as iREAD for identifying intron-retention similar to B above.

Is there a way, say samtools or other tool where I can find intron-exon or exon-intron read overhang extension.

As an example, image is given for easy understanding. Thanks for your help.

-A

Image showing intron-exon read overhang for top sample, whereas for the rest of samples below no such read overhang or coverage is observed.

![enter image description here][1]

[1]: http://funkyimg.com/i/2JmgG.jpg

Answer: Steps via `awk` and [BEDOPS][1].

Build a BED file of exons:

$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gtf.gz \
| gunzip -c - \
| awk '($3=="exon")' - \
| gtf2bed - \
| cut -f1-6 - \
> gencode.v28.annotation.exons.c1t6.bed

Convert transcripts to merged exons:

$ awk -f transcripts2mergedExons.awk gencode.v28.annotation.exons.c1t6.bed > gencode.v28.annotation.mergedExons.c1t6.bed

Make an exon-intron list:

$ awk -f mergedExons2exonIntronList.awk gencode.v28.annotation.mergedExons.c1t6.bed > gencode.v28.annotation.exonsAndIntrons.c1t6.bed

Convert these to junctions:

$ awk -f exonIntronList2JunctionList.awk gencode.v28.annotation.exonsAndIntrons.c1t6.bed > gencode.v28.annotation.exonsIntronJunctions.c1t6.bed

Optionally, pad them, e.g. by 25 nt around the junction (adjust as needed):

$ bedops --everything --range 25 gencode.v28.annotation.exonsIntronJunctions.c1t6.bed > gencode.v28.annotation.exonsIntronJunctions.pad25.c1t6.bed

Map reads to the padded junctions:

$ bedmap --echo --count --delim '\t' gencode.v28.annotation.exonsIntronJunctions.pad25.c1t6.bed <(bam2bed < reads.bam) > answer.bed

The file `answer.bed` will contain the junction and the number of reads that map to — overlap with — the optionally-padded junction, by one or more bases.

If you want the actual reads that map to the exon-intron junctions:

$ bedmap --echo-map --delim '\t' gencode.v28.annotation.exonsIntronJunctions.pad25.c1t6.bed <(bam2bed < reads.bam) | sort-bed - > answer.bed

In this second example, the file `answer.bed` will contain the reads that map to the junction, by one or more bases.

Some links to Github Gists with listed `awk` scripts:

* (`transcripts2mergedExons.awk`) https://gist.github.com/alexpreynolds/7f68953db11c7c00569141b577c01c81
* (`mergedExons2exonIntronList.awk`) https://gist.github.com/alexpreynolds/431216b2c3381b476bf91b149d35b2cc
* (`exonIntronList2JunctionList.awk`) https://gist.github.com/alexpreynolds/72c8deb54e85e4028698563d8c49f470

[1]: https://github.com/bedops/bedops
Question: In a previous question ["Code golf: mean length of fasta sequences"][1], Eric asked for some solutions `to get the average length of the sequences in a fasta file`.

I tried to anwser this question using the following [Erlang][2] code:

-module(golf).
-export([test/0]).

line([],{Sequences,Total}) -> {Sequences,Total};
line(">" ++ Rest,{Sequences,Total}) -> {Sequences+1,Total};
line(L,{Sequences,Total}) -> {Sequences,Total+string:len(string:strip(L))}.

scanLines(S,Sequences,Total)->
case io:get_line(S,'') of
eof -> {Sequences,Total};
{error,_} ->{Sequences,Total};
Line -> {S2,T2}=line(Line,{Sequences,Total}), scanLines(S,S2,T2)
end .

test()->
{Sequences,Total}=scanLines(standard_io,0,0),
io:format("~p\n",[Total/(1.0*Sequences)]),
halt().

**Compilation/Execution**:

erlc golf.erl
erl -noshell -s golf test < sequence.fasta
563.16

this code seems to work fine for a small fasta file but it takes hours to parse **uniprot_sprot.fasta** (in fact , I pressed Ctr-C). Why ? I'm an Erlang newbie, can you improve this code ?

[1]: /p/1758
[2]: http://www.erlang.org

Answer: The answer is [here](http://stackoverflow.com/questions/3296855/average-length-of-the-sequences-in-a-fasta-file-can-you-improve-this-erlang-co/3304288#3304288)
Question: Hello,

I have normalized gene expression data for 1095 patients with breast cancer
part of the data as follows

patient 1 patient 2 patient3 patient4 patient5 patient6

AASS	80.8588	135.3218	158.7152	20.7441	187.836	126.2016
AATF	2012.3344	990.1661	727.4445	1498.5344	1329.6371
AATK	179.534	209.8278	35.4275	13.5558	99.1263	51.2694
ABAT	2086.3408	77.9285	600.3779	101.0147	1564.1801 1439.816
ABCA11P	79.5614	91.9899	152.9806	65.3844	36.7641	85.9551
ABCA12	43.8556	1.0531	37.317	10.823	82.3253	4.9298
ABCA13	21.9278	2.6327	0.9447	0.9019	0.672	0

can I using this data find differentially expressed genes, what package in R can I used to get the DE genes?

Answer: If all of your samples have the same phenotype (cancer for example), why would you want to find DEGs?

You should have at least one control sample to compare with.

Or you want to see the overall expression pattern of your data?

In that case, you are going to need a lot of computer power.

Data with ~1,000 samples x ~30,000 genes is too big to run on an ordinary PC.
Question: What is the way to retrieve genomes from ncbi via biopython?
I am able to get a record for my genome of interest, and also i am able to download it manually from [search][1]

But how to download it inside a script?

from Bio import Entrez
Entrez.email = "my_email@email.ru"
handle = Entrez.esearch(db="genome", term="Drosophila eugracilis[Orgn]", idtype="acc")
record = Entrez.read(handle)
for i in record.keys():
print i,record[i]

[1]: https://www.ncbi.nlm.nih.gov/genome/?term=Drosophila%20eugracilis

Answer: `esearch/elink` commands would be different depending on the genome of interest.

For example, following command can be used to download assembly in the fasta format using NCBI Unix utils:

esearch -db genome -query "6863"|elink -target assembly|esummary|xtract -pattern FtpPath_RefSeq -element FtpPath_RefSeq

This should give you the URL where assembly files are saved and relevant files can be downloaded from that URL.
Question: Hi,

I am working in R environment and I have a long list of gene names which I am showing few first objects of it here:

```
[1] SPNCRNA.1436,omh5,snR95
[2] snR46
[3] snR10
[4] SPNCRNA.1651,SPNCRNA.515
[5] snR42
[6] SPNCRNA.1094,SPNCRNA.1095,SPRRNA.47,SPRRNA.48
[7] snR88
[8] SPNCRNA.497
[9] SPSNORNA.54
[10] snoR39b
```

I am wondering if there is any way to split the indexes with several gene names in them into individual ones? The list is so long.

Thanks!

Answer: You can `sapply` with `strsplit` to split each element of the vector by `,`

Edit: Sorry, I was half asleep when I wrote this answer, so could not test code before posting it here. I've now added the logic behind my solution, code and output.

Logic:

Apply `strsplit` on each element of the vector and flatten the resulting list of vectors using `unlist`.

Code:

```
> x <- c("SPNCRNA.1436,omh5,snR95","snR46", "snR10", "SPNCRNA.1651,SPNCRNA.515","SPNCRNA.1094,SPNCRNA.1095,SPRRNA.47,SPRRNA.48")
> listOfSplitStringVectors<-sapply(x,function(i) strsplit(i,",")
> flattenedVectorOfNames=unlist(listOfSplitStringVectors,recursive = TRUE,use.names = FALSE)
```

Output:

```
[1] "SPNCRNA.1436" "omh5" "snR95" "snR46" "snR10"
[6] "SPNCRNA.1651" "SPNCRNA.515" "SPNCRNA.1094" "SPNCRNA.1095" "SPRRNA.47"
[11] "SPRRNA.48"
```
Question: Hello!
I am trying to use RNA-seq data to accomplish two objectives in a case vs. control study:

1 - Verify if the expression level is distinct between the case sample and the controls.

2 - We need to identify DE genes between conditions.

Unfortunately, we have only one case with the mutation, but we have six biological replicates for the control group. Is it ok to use the package edgeR, or we need to try something else?

Thanks!

Answer: There's no problem with this. edgeR handles unequal sample sizes automatically. It will simply estimate variability (i.e., the dispersion parameter) from the controls, where you do have replicates, and will then apply the same dispersion values to both cases and controls when conducting tests and computing fold changes.

This is a frequently asked question:

- https://www.biostars.org/p/266211/
- https://support.bioconductor.org/p/113521/
- https://support.bioconductor.org/p/76631/
- https://support.bioconductor.org/p/63585/
- https://support.bioconductor.org/p/61904/
- https://support.bioconductor.org/p/53627/
Question: Hi All

I am confused in extracting the RGID, RGLB and RGPU from this read, any help, thanks?

@E00461:116:AU170602261:1:1101:20456:1415 1:N:0:NTCCAGGT

And if I have each sample is a tumor from a different patient, should I use RGSM as the different sample patient ID?

Thanks

Answer: Hi Tania

I used this tutorial to get those values from the fastq files, https://angus.readthedocs.io/en/2017/GATK_pipeline.html

Jump to: **4. Add Read group information and do mapping**
Question: is it a good practice to remove all the unmapped reads with:

samtools view -F 4

after they've been mapped with bwa ? The bam files would be smaller and the remaining operations would be faster isn't it ?

or shall I regret it later ?

Answer: To future-proof your data, it seems reasonable to hold on to the unmapped reads. For one, future algorithms may do a better job and allow you to recover some of that data. For another, future genome builds may resolve poorly assembled regions and allow additional reads to be mapped.

Neither of these improvements is likely to enable huge discoveries, but the cost you're paying in storage is pretty minimal, compared to the costs of sample collection and sequencing. The speed hit probably isn't as bad as you think either, since the bam is indexed. Smart algorithms will make use of that information and not even have to consider those unmapped reads.
Question: I'm working on a dataset in an insect species. I was using [KAAS][1] to annotate the genes. In the parameter for "GENES data set", I only selected some species from insects.

But in the results, the Pathway Mapping gave pathways definitely human specific, E.g. "Coronavirus disease - COVID-19". In this pathway, there is one KO [K02580][2]. This KO belongs to many pathways. One of them is "Coronavirus disease - COVID-19". My question is:

Is there a way that I can exclude some of the pathways that we know definitely not related to insects? Any suggestions? Thanks.

[1]: https://www.genome.jp/kaas-bin/kaas_main
[2]: https://www.genome.jp/dbget-bin/www_bget?ko:K02580

Answer: Yes I do. This is what I did:
1. In the webpage of KAAS, I'll get a list of species I would like to limit the pathways to. eg. dme, dya.
2. Download the pathways for the species you selected using the link below:

[http://rest.kegg.jp/link/pathway/dme][1]

[http://rest.kegg.jp/link/pathway/dya][2]

In the second column, you will see the IDs. I put an example down below here:
dme:Dmel_CG10160	path:dme00010

The idea is that dme00010 and ko00010 are referencing to the same pathway:

[https://www.genome.jp/entry/dme00010][3]

[https://www.genome.jp/entry/ko00010][4]

So you can only include that pathways that are in the selected species based on the IDs without letters prefixed.

Hope this helps.

[1]: http://rest.kegg.jp/link/pathway/dme
[2]: http://rest.kegg.jp/link/pathway/dya
[3]: https://www.genome.jp/entry/dme00010
[4]: https://www.genome.jp/entry/ko00010
Question: Hi all,

I have paired-end 4 ATAC-seq data (2 replicates for 2 samples). I have done aligning using Bowtie2. I did filter MT reads and duplicates using Picard, then performed peak calling on Bam file using MACS2. Also I did differential peak analysis using deeptools and filter them by FDR<0.05 and abs(2foldchange)>2.

After these, I generated density peak heatmaps using deeptools. However, on the figure top on the heatmaps height of peaks are not the same for 4 files although I normalized bam files while converting to bigwig using bamCoverage.

My questions are:
Should the height of those peaks be the same or slight change is acceptable?
If not how can I normalize the data? Should I normalize bam files then do the peak calling again if so which tool you suggest? or diffbind normalization okay?
Also, I am really confused about the coverage file normalization and peak normalization.
Lastly, as written in this post https://www.biostars.org/p/308976/ how can I downsample each sample?

If you could explain these I will appreciate it.

![peak signal heatmap][1]

[1]: https://www.dropbox.com/home?preview=image.png

Answer: In some cases normalizing only for sequencing depth might be enough. Often it is not due to differences in library composition and different signal-to-noise ratios. I prefer to scale my bigwig files (or whatever counts you want to normalize) with the normalization factors from `edgeR`. Code examples and some details in the linked thread: https://www.biostars.org/p/413626/#414440
Question: So I have a sam file with extended CIGAR format like this:

17H87=1D12=1D2=3D32=1D4=2D5=2D13=....

And I want to convert it to a sam file containing regular CIGAR format that will look like this:

17S83M1D14M1D4M3D31M1D2M2D6M...

(they do not represent the same reads... I want only to visualize the problem)

Anyone knows some software/script that will convert sam/bam with extended CIGAR to a sam/bam with regular CIGAR?
It seems like a standard problem that might already have a solution so I want to ask here before solving it the hard way(coding)

Answer: `reformat.sh` from [BBMap suite][1] can do this: `reformat.sh in=orig.sam out=new.sam sam=1.3`

[1]: https://sourceforge.net/projects/bbmap/
Question: Hi all,
I've just written a tool adding one or more **extra column** in a [VCF file](http://www.1000genomes.org/wiki/doku.php?id=1000_genomes:analysis:vcf3.3). The header now looks like this:

(...)
#CHROM POS ID REF ALT QUAL FILTER INFO MY_COL1 MY_COL2 FORMAT NA00001 NA00002 NA00003
(...)

Is there something in the **VCF** spec saying that another column can't be added ?
because when I used **VCFTOOLS**, it says:

vcftools --vcf file.vcf
(...)
Scanning file.vcf ...
Ninth Header entry should be FORMAT: MY_COL1
Currently scanning CHROM: 19
Currently scanning CHROM: 20
Currently scanning CHROM: X

Answer: BEDTools now supports VCF and you can tack on any number of columns you want. That said, if you are looking for specific functionality within VCFTools, then this isn't helpful at all.
Question: Hello,

I'm trying to build my first Snakemake worklflow. I've got four sets of short reads for two samples (left/right reads for samples HB001_naive/coevolved) that I'd like to correct with a script bbduk.sh. I'm trying to write a rule for that:

```
SAMPLES=["HB001_naive", "HB001_coevolved"]
READ_SETS=["1", "2"]

rule all:
input:
expand("1_reads_qc/{sample}/{sample}_{read_set}.bbduk.fq.gz", sample=SAMPLES, read_set=READ_SETS)

rule bbduk:
input:
expand("0_reads_raw/{sample}_{read_set}.fq.gz", sample=SAMPLES, read_set=READ_SETS),
output:
expand("1_reads_qc/{sample}/{sample}_{read_set}.bbduk.fq.gz", sample=SAMPLES, read_set=READ_SETS),
conda:
"envs/bbtools.yaml"
threads: 12
shell:
"bbduk.sh -Xmx60g t={threads} \
in1={input} in2={input} \
out1={output} out2={output} \
ref=/home/scro4331/.conda/envs/bbtools/opt/bbmap-38.18/resources/adapters.fa \
ktrim=r k=23 mink=11 hdist=2 maq=10 minlen=100 tpe tbo \
stats=bbduk.contaminants"
```
This lists all the input files after `in1=` and `in2=`, same for output. However, I'd obviously like Snakemake to consider first only the naive read sets (as `in1=0_reads_raw/HB001_naive_1.fq.gz` and `in2=in1=0_reads_raw/HB001_naive_2.fq.gz`), then the coevolved ones.

Any suggestions on how to make it work, allowing flexibility in adding more samples? I'd appreciate any help!

Answer: By doing

expand("0_reads_raw/{sample}_{read_set}.fq.gz", sample=SAMPLES, read_set=READ_SETS)

You're passing at the same time every permutation of `$sample` and `$read_set` , creating a list for the variable `input`. Doing that inside your rule you're passing this file list two times, for both `in1` and `in2`.

Better splitting input/output into the two pair

rule bbduk:
input:
i1 = "0_reads_raw/{sample}_1.fq.gz",
i2 = "0_reads_raw/{sample}_2.fq.gz"
output:
i1 = "1_reads_qc/{sample}/{sample}_1.bbduk.fq.gz",
i2 = "1_reads_qc/{sample}/{sample}_2.bbduk.fq.gz"
shell:
"""
bbduk.sh -Xmx60g t={threads} \
in1={input.i1} in2={input.i2} \
out1={output.i1} out2={output.i2} \
ref=/home/scro4331/.conda/envs/bbtools/opt/bbmap-38.18/resources/adapters.fa \
ktrim=r k=23 mink=11 hdist=2 maq=10 minlen=100 tpe tbo \
stats=bbduk.contaminants
"""
Question: Hi all

I am using cytoscape to evaluate the overlaps between multiple PPi networks.

After importing the networks from a txt file I used the "merge" tool using the "union" option to create the merged network. This part works well but ideally I would like to color differently each network and highlight the overlapping nodes with different colors/shapes, but I can't seem to find the way to do it. I tried to google the various "visualize/color overlapping nodes", someone suggested a plugin [here][1] but I can't find any indications that this was actually done.

Any suggestions/workarounds?

Thanks :)

[1]: http://wiki.cytoscape.org/Future_Cytoscape_Features

Answer: Create a mapping file in which you assign to each node the category it belongs to (i.e. network 1, network 2 or shared). Upload this file into cytoscape as a table using node names/IDs as keys. This will make another column in the node table associating the category as an extra node attribute. Then in the style tab, assign a color to each category using the discrete color mapping option of the node color.
Question: I concatenated two vcf files by scripting, keeping only the header from the first file and adding an extra column in the end to know which file each line is coming from.

I do this:

tabix -p vcf f1PotMut.vcf.gz

and get this:

[ti_index_core] the chromosome blocks not continuous at line 87352, is the file sorted? [pos 131430]

so I do this

vcf-sort f1PotMut.vcf.gz

but the problem remains. I get the same error when trying to tabex my "sorted" file. I somehow feel the vcf-sort command did not work. Any help?

Answer: above you are not showing but for the sake of covering that: are you redirecting the output of the vcf-sort into a different file then are you running tabix on this new file?
Question: I have some Genomic coordinates based on hg19 which is equal to GRCH37 and I would like to analyze these coordinates using actual genome sequences. In Ensembl, There are several releases of Grch like GRCH37.70 and GRCH37.75 etc. I am wondering what is the actual difference between all these releases of GRCH37 ?! Are the coordinates the same for all these releases?!

Answer: This is the website of the Genome Reference Consortium, might be informative for you to poke around there

hg19 is typically thought of as a subset of GRCh37 with the GRCh37p# representing different minor / patch releases (changes do not affect coordinates)

Ensembl uses a different versioning system, you can find their list of changes here: http://www.ensembl.org/Help/ArchiveList though it seems like the latest release is mostly schema changes with some changes to annotations.
Question: Hey guys:
I am doing **RNA-seq** analysis and it seems that the quality of my reads is not desirable.
Below is **a typical fastqc report** for my data.
I have read many tutorial about fastqc, from my understanding, **it seems that the 1-10 bp are adaptor sequences**. But in the adaptor content section, there is no waining.

----------

- I am wandering if my understanding is right?
- Should I use **trimmomatic** to cut adaptor sequences?

![Failed per base sequence content][1]

[1]: https://lh3.googleusercontent.com/Xp-_0jmKsw8p2MVrOx_X-hJbHIiKBfUFZ700lNT5nv0MKlma6DSgYonxLAbkHzkYNEaWHszNcvI6RjCPhk6YghTHfI7dT1_2GUyzPr4I5MdEnee1_qTSB46nbvwTJFh_TEASRgS1MEGIawbdc7zZpnWWd_bFL_CeNB0-MBk6PePyflSrMjeLZFlFu-1onTqcypuSMKeYw5Wty0c-pgPSCHWIHkTxrg3d0Y-XUhztxFRETL2sFQQXjtSGmZ5PhguduitzrZcnOijNlVK5KK_0hZ5Aqd-3pf7OHjJbRXfP02dxda5HQOrALIo-_DYoHszlbMDSFeLwSQerz80bPd2CNEsLmE79KsYBSqjOIeg_aH8npy70IlZc2pT34_pZul8pZapn_Yb_QN2tiazcaPPeYJvC-uVttpkoepsoaXT7Oz5lbB1OBkCrx1Ye3dyzpAx_bGhvdnnIXhOFYKbI44Nr1mh1ns7up0uj6uukZDqFOvRAkEFRNeNPHuqZlxGFfnCY1qAuM3Kqy3c0UjhVdDmFlTNVnR6wICnRR8v4EL7HxCbdPAnZ7JlHHzpaPHtKTxlMG80d51jSwWVLPa1GnaveRBsjfp3Ek5NVZ8bhv5UxhQXhPKRi8L3XBZ6-bg0t3MGomxWRSaE54EN6HD2GV1NKXE8muzTGcw=w1407-h809-no

Answer: Hi,

This is totally normal for RNA-sequencing data even after removal of the adapter sequences.
The random hexamer primers, which are used to generate the cDNA library from your RNA transcripts were shown to not bind completly random. This non-random binding leads to this bias in "per base sequence content" from base 1-15.

You don't have to trim these sequences, if fastqc does not report recognizes adapters. Just be aware, that the used primers do not lead to a completly random amplification.
Question: Could I have specific row labels in bold text but not all of them when creating a heatmap using R? (eg. using heatmap.2 function)

Answer: I just figured out that you can also give labels different colors using *heatmap.2* with argument ***colRow***.
Question: Hi Biostars,

I have a set of about 6,000 groups of orthologous proteins. Each orthologous group has a representative protein from anywhere between 7 and 16 species (1 per species). I'm trying to do ancestral state reconstructions for each of these groups based on their well-established phylogenies, and then to compare each amino acid position in a species of interest to the reconstructed sequences at key ancestral nodes.

I've been trying to use PAML (codeml) for this purpose, following the method used in this blog post http://evosite3d.blogspot.com.au/2014/09/tutorial-on-ancestral-sequence.html

My problem is that when using clock = 0 (no molecular clock), PAML requires an unrooted tree. According to PAML's manual "a rooted tree has a bifurcation at the root, while an unrooted tree has a trifurcation or multifurcation at the root."

This is problem is that the mammalian tree that I have (composed of eutherians and marsupials) is a well-established bifurcating tree. How can I do an ancestral state reconstruction with PAML which requires a multifurcating tree when the true phylogeny is bifurcating? Adding an outgroup like platypus or chicken would still be a rooted tree since they're both outgroups and a polytomy of chicken/platypus, eutherians and marsupials would be false. I'm sure I've just deeply misunderstood something along the way. Any help would be really greatly appreciated!

Answer: The short answer is to unroot the tree. This can be accomplished easily in R using `unroot()` in `ape`.

It is possible for a tree that is normally bifurcating when rooted to have a node that has an order greater than two when unrooted. Consider a three-taxa case. There is one possible unrooted tree (shaped like a 'Y') and three possible rooted trees (whether rooted using the first, second, or third taxon) that are bifurcating. Remember that the placement of the root is a hypothesis.