cannin/biostars_qa huggingface dataset:
Question: Hi guys,

I need to identify insertions of specific sequences in my genomes. I have WGS data with read length of 250 nucl. Previously I was working on identification of large sequences insertions and simply looked for discordantly mapped reads for that using bowtie2 and bbmap.
Now I am looking for sequences that are shorter than read length. As I want these data to be comparable to my previous computations I want to use bowtie2 for alignment again. As far as I understand, I should use local mode. I've tried this, but got alignment of very short parts of my reads. Is there a possibility to require the alignment length to be equal to reference sequence length? When I identify the reads that are mapped to the reference sequences, I will find where are their pairs mapped to the reference genome to find the insertion site, is it a correct strategy?

Thanks id advance!

Answer: I'd recommend using lumpy or delly2 for small indels, rather than trying reinvent your own method

https://github.com/dellytools/delly

https://github.com/arq5x/lumpy-sv
Question: I am new to plink and have a (hopefully) simple data management operation, but little idea of how to implement it. The essentials are listed here:

1. I have three genotype call files: `data.bed`, `data.fam`, and `data.bim`.
2. I created a file `snps.txt` containing a list of markers.
3. My goal is to create a file with subject IDs listed in rows and SNP IDs (as per `snps.txt`) as column headers with genotypes listed underneath.

Question: How might I go about performing this task?

Answer: If you're fine with the main body of the file containing allele counts (0/1/2) rather than allele names, you can use

plink --bfile data --extract snps.txt --recodeA

and then use Unix cut to remove any header columns you don't want.

Otherwise,

plink --bfile data --extract snps.txt --recode compound-genotypes

(this requires PLINK 1.9) will almost get you there, but you'll need to add a header line on your own. Or

plink --bfile data --extract snps.txt --recode

works if you want two columns per SNP.
Question: I am attempting to run a comparison a batch of alignments (bam format) to assess spread. However, some of my datasets are larger than others and I wish to subsample down to an equal number of reads per sample.

If I wanted to randomly extract 1 Million reads from a bam file is there a method to do this?

Note: I am fully aware of the samtools and picard methods which allow you to reduce by a proportion (i.e. the flag `-s 0.33`) but that would not result in an fixed number of reads per sample, which is what I need, but a reduced proportion *per sample* which doesn't help.

Bam subsampling has previously been talked about it but always the proportional data reduction, not the fixed number required: https://www.biostars.org/p/44527/, https://www.biostars.org/p/76791/#76796, https://www.biostars.org/p/110041, https://broadinstitute.github.io/picard/command-line-overview.html

Edit: Also, I've come across [bamtools][1] but haven't been able to get it to work and seems to have not been updated in quite some time

[1]: https://hcc-docs.unl.edu/display/HCCDOC/Running+BamTools+Commands

Answer: This is based on some of the earlier threads mentioned in the comments, but you can use the GNU `shuf` utility to get a fixed number of random lines in a file. [Recent versions][1] of the tool use [reservoir sampling][2], which will improve efficiency. You could run something like:

shuf -n 1000000 input.sam > output.sam

To get the header correctly, you could do something like:

cat <(samtools view -SH input.sam) <(samtools view -S input.sam | shuf -n 1000000) > output.sam

[1]: http://lists.gnu.org/archive/html/coreutils/2013-12/msg00166.html
[2]: http://en.wikipedia.org/wiki/Reservoir_sampling
Question: Hi,

I've a small question concerning differential expression analysis. So I've a 12 samples divided into three groups (A, B and C ). If a do that :

samples <- read.table("samples.txt",sep="\t",as.is=T)
colnames(samples)<-c("sampleName","fileName","condition")
dds1 <- DESeqDataSetFromHTSeqCount(sampleTable=samples,directory="./",design= ~ condition)
dds1 <- DESeq(dds1)
res <- results(dds1)

What will the p-value represent? If it's significatif, is it because the gene of interest is up or down regulated in one of the three groups?

To do pairwise comparison, I do that. Is it correct?

res.A.B <- results(dds1, contrast=c("condition","A","B"))
res.A.C <- results(dds1, contrast=c("condition","A","C"))
res.B.C <- results(dds1, contrast=c("condition","B","C"))

Thanks

Answer: By default, DESeq2 uses a Wald test to compute the p-values, so it's likely that `res` just has the output of the A vs. C contrast (rather than using a likelihood ratio test to compare `~condition` versus `~1`). Regardless of the test, the direction of change and the p-value are separate, so you'll have significant findings in both directions.

Yes, the method you showed for the pairwise comparisons looks correct.

BTW, if you really do want to test whether something is changed due to condition in general (i.e., do something like a one-way ANOVA), you can just use the LRT methods. Keep in mind that the direction of change there can be questionable, since you're then asking if the full model fits better than the reduced rather than asking of parameter estimates (i.e., fold-changes) are different than 0 (or a different threshold, should you specify something else...though that's not common).
Question: Hi All,

I need to do some permutation test to random sampling millions of SNPs set genotyping data from 1000 genome phase III dataset. each SNP-set only contains ~300 SNPs, however, I found the extract process is quite slow. for example, for the ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf, to extract these 300 SNPs require almost 20 mins.

plink --vcf ~/hpc/db/hg19/1000Genome/chr6.uni.vcf --extract chr6.rs6457620.txt --r inter-chr dprime --out rs6457620

Is there any fast way to extract the subset quickly?

real 17m52.822s
user 8m39.925s
sys 1m20.681s

Done! Thanks.

```
tabix ~/hpc/db/hg19/1000Genome/chr6.uni.vcf.gz -R chr6.rs6457620.txt > output.vcf
```

100 times faster!!!

Thanks.

Answer: Sounds like a job for [Tabix][1]. Have a look at the manual and the `-R` option. It allows quick retrieval of regions from a tabix-indexed VCF, GFF/GTF/BED file.

[1]: http://www.htslib.org/doc/tabix.html
Question: I've been using DAVID for gene ontology IDs for years but this -> https://www.biostars.org/p/245/ says it hasn't been updated in like 6 years. Is there anything else that's useful and up-to-date? I usually use Uniprot, KEGG, and GO IDs for everything. http://geneontology.org/ wants me to input a species but the species I deal with aren't really model organisms. I have a lot of random microbes.

**What can I use instead of DAVID? Preferably accessible through Python or web interface**. If there is some REALLY good that's for R I'm open but that's a last resort for me (It takes me a while to do anything in R).

Answer: The Gene Ontology consortium has created their own enrichment tool through AmiGO 2 http://amigo.geneontology.org/amigo (and/or with Panther http://pantherdb.org/webservices/go/overrep.jsp). This will be *the* most up-to-date tool you can use.

Intermine provides up-to-date enrichment tools (web interfaces or accessible by API) for a lot of species here: http://intermine.org/
Choose the closest species mine if you can, it doesn't need to be a perfect match - especially for GO which is species-neutral; though some of the other functions they index can (or should) be species-aware. For microbes, YeastMine/Saccharomyces cerevisiae database might be the best? InterMine/MODs get updated GO annotations as they come out, or very soon after (~ up to 3 month delay depending on production cycles).
Question: According to the spec a template is "A DNA/RNA sequence part of which is sequenced on a sequencing machine or assembled from raw sequences." But what does that mean? The first field "QNAME" is explained as "Query template NAME." I know that usually means read name, so what does "template" have to do with it? Also from the spec, "RNEXT: Reference sequence name of the primary alignment of the NEXT read in the template. For the last read, the next read is the first read in the template. " Does template just me "file" then? But then a template can be "single-segment," in which case TLEN is set to 0. Try as I might, I can't make sense of that. Your explanation will be much appreciated.

Answer: In most cases, templates are reads or pairs of them. However, you could imagine a case where you assembled reads into contigs and then aligned those against a reference genome. Then the template would be the contig.

RNEXT is the name of the chromosome or contig to which the next template in a pair (or larger group) aligns.
Question: When using PRSice for PRS calculation of target data using a file of variants from the PGS catalog, does the variant file from the PGS catalog replace the GWAS summary statistics (referred to as "base data" in the tutorial)?

I had previously found this similar question: https://www.biostars.org/p/407615/, but just making sure that the list of PRS variants replaces the summary statistics/base data.

Thanks in advance.

Answer: As long as it has the SNP, Effective allele, effect size and p-value column, then it is fine. If there isn't a p-value column, you can sometimes make a "fake" p-value column with all 1s or all 0s. Depending on the data, you might also want to skip clumping and thresholding.

Sam
Question: Hi,

Does anyone has detailed documentation of PRSice-2 software? The official website has been down for several weeks.

Thank you!

Mengna

Answer: We lost the domain name because we forgot the login password... You can now revert to the repo doc (which is the same website as before, just that previously it was redirected to PRSice.info)

Sam
Question: siRNA's are 19-21 long RNA duplexes with 2 base 3' overhangs present in the RNA-interference pathway.
I would like to calculate the local internal stability of several RNA duplexes in siRNA, as described in Khvorova et al. (2003) based on (delta G) calculations according to Freier et al (1986). I also want to calculate local stability at other temperatures than 37 degrees (not every bug lives at 37 degr. C) to achieve plots similar to Figure 2 in Khvorova et al. How well does the stability calculation transfer?

In Kvorova et al., they used the proprietary OLIGO software for this purpose. Is there an open-source alternative to this procedure (Khvorova et al.):

The Oligo 5.0 Primer Analysis Software (National Biosciences, Inc., Plymouth, MN) was used for calculations of internal stability profiles using the free energy values in Table 1 (−kcal/mol, adjusted to 25°C).

... The internal stability values reflect the stability of pentamer subsequences within the sequence under investigation and were calculated according to the nearest neighbor method (Freier et al., 1986). For the calculation of the average internal stability values of the terminal four bases on the 3′ end of the molecule, the 19 bp targeting sites were extended based on the target mRNA or miRNA precursor sequences.

I guess the calculation should be rather basic, therefore a formula for the nearest neighbor method of Freier should be sufficient.

Sorry if this is a very basic and partially biochemistry/biophysics question.
Edit: Seemingly, it is not at all as easy as I thought, even though any primer designer software and RNA-folding software should internally compute these stabilities. I tried for example Pairfold and mfold, but that outputs only minimum free energy (which is 0 for perfect pairing). Another hint: the claculation should be similar as for "3' Max stability" in eprimer3 but with thermodynamic parameters adjusted for RNA.

Answer: Have you tried RNAfold in Vienna RNA?

main page: http://www.tbi.univie.ac.at/~ivo/RNA/

latest version: http://www.tbi.univie.ac.at/~ronny/programs/vrna2g.html

-T allows temperature scaling.
Question: I constantly download tools and programs to work with. I don't want to have to create a new PATH every time. How can I organize ALL my tools/ programs/ scripts (old and new) into one place so that I don't have to worry about running commands from anywhere on my terminal? any ideas?

Answer: Personally, I try to avoid `make install` from any software like the plague. Usually, I forget about the software so it just ends up filling my disk space and leading to other issues (not receiving any updates etc.). Instead these days I have a directory with a bunch of 3rd party software source files (commonly also git repo clones). There I build the tools and have a symbolic link from `$HOME/bin` to each tool with `$HOME/bin` being in my `$PATH`. That way uninstalling is easy, updating is easy (`git pull && make`), and I even have the source code at hand if I stumble upon a bug. This scheme is almost as good as making my own package via [makepkg](https://wiki.archlinux.org/index.php/Makepkg) or smthn similar.
Question: I want to plot the number of positions in a sliding window of **1000** and a step of **20** for each sample **(A-D)**.

Interpretation:
* `1`: position exists;
* `NA`: position does not exist.

I have tested a dozen tools in `bash`, `R` and other but I am looking for a `Python` solution.
Your advice please.

#This is an example of my data:
window = 1000
step = 20

# Example of dataframe
POSITION	A	B	C	D
1250	1	1	1	1
1750	NA	1	NA	1
1786	1	NA	1	1
1812	1	1	1	1
1855	1	1	1	1
1896	1	NA	1	NA
2635	NA	1	1	1
1689	1	1	NA	NA
3250	1	1	1	1
3655	1	NA	1	1
3589	NA	1	1	1

I am looking for some thing like this:

![Polymorphism density plot using a window size of 1,000,000 with an increment (step) of 100,000][1]

Any help will be appreciated!

[1]: /media/images/a29910c4-0547-4588-9d66-48de0400

Answer: This is a Python solution:

import pandas as pd
import numpy as np
import seaborn as sns
df = pd.DataFrame({'POSITION': [1250,
1750,
1786,
1812,
1855,
1896,
2635,
1689,
3250,
3655,
3589],
'A': [1.0, np.nan, 1.0, 1.0, 1.0, 1.0, np.nan, 1.0, 1.0, 1.0, np.nan],
'B': [1.0, 1.0, np.nan, 1.0, 1.0, np.nan, 1.0, 1.0, 1.0, np.nan, 1.0],
'C': [1.0, np.nan, 1.0, 1.0, 1.0, 1.0, 1.0, np.nan, 1.0, 1.0, 1.0],
'D': [1.0, 1.0, 1.0, 1.0, 1.0, np.nan, 1.0, np.nan, 1.0, 1.0, 1.0]})

window = 5
step = 2

df = df.set_index('POSITION').rolling(window).count().reset_index().iloc[::step, :]
df = df.melt(id_vars='POSITION', value_vars=['A','B','C','D'], value_name='polym', var_name='chromop')
sns.lineplot(data=df, x='POSITION',y='polym',hue='chromop')
Question: Hi there,

Fairly new to this area so will try to explain as clearly as possible. I have two RNA-Seq data sets - one corresponds to a series of cancer cell lines, the other to cell lines we are using as a model of 'normal' epithelium. The expression units of one data set are in FPKM the other RPKM. I think I superficially understand the difference between these units, in that one is used in mapping transcripts in single end sequencing the other paired end sequencing.

My question is are these units directly comparable? The analysis I wish to carry out is fairly straight forward - I have a predefined list of around 500 genes, and simply want to compare differences in expression between the non-cancer/cancer background, in terms of which of these genes they are expressing at all as well as the relative expression levels. For which transcripts are expressed I had intended to use any value over 0 (FPKM or RPKM) as denoting expression of a transcript, but am unsure if I can compare relative abundance.

I should add that I only have access to the raw data of one dataset, the other is as a results table sent by collaborators.

Answer: In theory RPKM should just be twice FPKM, but that is unlikely to actually be the case, as is noted by the other commenters, differences in processing protocol, are likely to mean they are not comparable, and would not be even if they were both FPKM or both RPKM.

You **might** be able to get some meaning from the dataset if you are just doing presence/absence calls, but I wouldn't use a 0 threshold. Instead convert your values to TPM (transcripts per million) but normalising each value to the sum of expression for that sample and multiplying by 1 million. E.g. if I had RPKMs from 4 genes of 0, 1, 5 and 4 the TPMs would be (0*x10^6)/10, (1*10^6)/10, (5*10^6)/10 and (4*10^6)/10. Transcripts Per Million is much easier to reason with and more comparable between samples. As a very approximate rule of thumb, a human cell has around 200,000 mRNAs in it at any one time, so a TPM of 5 corresponds to one molecule per cell and can serve as a handy rule of thumb for calling present/absent.

Be careful though as differences in processing might lead to false positive differences - e.g. if one set was processed allowing multimapping and the other not, then you will be a present call for paralogs in one sample and an absent call in the other, irrespective of whether there is actually an expression difference.
Question: Hi
I am new in python and want to see how this can be done in python (I can do this in R). I have a text file `myfile.txt` with one column and thousands of rows as shown below. I want to convert this to fasta `result.fasta` format as shown below. How can I do this in python?

`myfile.txt`

ATGTGTGGTTTTCCCCC
ATTGGCGGGGTTTTTCAGGGG
ATGGGGGGGCCCCCCCCAAAAAA
TTGGTGGGGGGGGGGGGAA

`result.fasta`

>1
ATGTGTGGTTTTCCCCC
>2
ATTGGCGGGGTTTTTCAGGGG
>3
ATGGGGGGGCCCCCCCCAAAAAA
>4
TTGGTGGGGGGGGGGGGAA

Answer: Create a new script called `ConvertFASTA.py`:

import sys

#File input
fileInput = open(sys.argv[1], "r")

#File output
fileOutput = open(sys.argv[2], "w")

#Seq count
count = 1 ;

#Loop through each line in the input file
print "Converting to FASTA..."
for strLine in fileInput:

#Strip the endline character from each input line
strLine = strLine.rstrip("\n")

#Output the header
fileOutput.write(">" + str(count) + "\n")
fileOutput.write(strLine + "\n")

count = count + 1
print ("Done.")

#Close the input and output file
fileInput.close()
fileOutput.close()

Then, run it with:

python ConvertFASTA.py myfile.txt result.fasta
Question: Hi all,

I recently downloaded HMMER in order to use **hmmer search** to identify protein sequences.

I know its possible to do your own database in hmmer but I would like to use the Pfam database which is way more complete.

First, can the Pfam database be downloaded? And if so, where? I can't seem to find it online...
Then, how large is it?

Thank you so much in advance for your help!

Have a great day

Answer: PFAM files: ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release

Pfam-A.hmm - Pfam-A HMMs in an HMM library searchable with the hmmscan program.
Pfam-A.hmm.dat - Data associated with each HMM required for pfam_scan.pl
Question: Im working on some analysis with miRNA expression data and miRNA targets information from different sources like miRBase, Targetscan etc...

The names of miRNAs in the expression data and miRNAs targets information data need to be same for my analysis

But I see that in both of these data, names look different.

For eg: In miRNA expression data, there is a miRNA `hsa-miR-373` and in miRNA target information I see there are `hsa-miR-373-3p, hsa-miR-373-5p`.

similarly, In expression data `hsa-mir-101-1 and hsa-mir-101-2` and miRNA target information it is `hsa-mir-101-3p`

Question: Can I change the names of `hsa-miR-373-3p, hsa-miR-373-5p` into `hsa-miR-373`?

Answer: miRNA nomenclature can be a bit confusing. The direct answer to your question is yes somtimes, but not always, and not without checking.

miRNAs are expressed as hairpins. Thus, the precursor hairpin to hsa-miR-373 looks like:

- g uuuuug |
5' gggauacuc aaaau ggggcgcuuucc u |
||||.||.| ||||| |..||.|||.|| | - This is pre-hsa-miR-373
3' cccuguggg uuuua cuucgugaaggg c |
g g ucaugu |

This is processed to make a double stranded RNA, which contains two short RNAs that could potentially act as miRNAs:

acuc-aaaaugggggcgcuuucc <---- This is hsa-miR-373-5p
||.| ||||| |..||||||.
ugugggguuuuagcuucgugaag <---- This is hsa-miR-373-3p

For most RNAs one of these arms is unstable and is quickly degraded, so it doesn't get chance to have any activity. The other is stable and get incorporated into RISC - the enzyme that does the business for miRNAs. In the case of hsa-miR-373 the 3p arm is stable and the 5p arm is unstable. Thus the 3p arm is also known as hsa-miR-373 and the 5p arm is also known as hsa-miR-373*. You can't treat hsa-miR-373-5p/miR-373* as the same as miR-373, because it has an entirely different sequence, and would have different targets if it were ever stable enough to get incorporated in to RISC. The information on which arm is which is in mirbase. There are some miRNAs however, where both arms are stable and must be treated as independent miRNAs.

The next complication is that some miRNAs appear more than once in the genome. hsa-miR-101 is an example of this. hsa-miR101-1 is located on chromosome 1, and hsa-miR-101-2 is located on chromosome 9. But both encode miRNA precursors that have identical 3p arms. That is the sequence of the 3p mature miRNA product produced from both the hsa-miR-101-1 and hsa-miR-101-2 loci is `uacaguacugugauaacugaa`, and this is the stable arm. If you sequence this, you have no idea whether it came from miR-101-1 or miR-101-2. However, the unstable 5p arms are different for miR-101-1 and miR-101-2. So officially there are 3 mature miR-101 RNAs in humans - miR-101-1-5p, miR-101-2-5p and miR-101-3p, expressed from two loci. In this case you can safely treat miR-101-3p as miR-101 irrespective of which locus it came from. Sometimes both arms are identical. Again, the information on which mature products are the same between the two copies, and which is the stable and unstable arm, is on mirbase.