cannin/biostars_qa huggingface dataset:
Question: Hi,

My question may seem so simple. Could you tell me what is the difference between aligning and mapping the short reads to the reference genome?

And also what is the difference between Pairwise alignment, Multiple sequence alignment and Short-Read Sequence Alignment?

Thanks

Answer: For your second question, pairwise alignment (e.g., Smith and Waterman) is between two sequences and multiple sequence alignment between more than 2 sequences (e.g., clustalW). Short read aligners are usually pairwise!
Question: This seems simple maths behind, But I always get confused with the terminology of **"Fold change", "Log Fold Change" *(LFC)*, and "Percentage of Change"**., When I see some of the figure descriptions given in **volcano plots**, specially X axis, I get confused.

i.e. [Page 64][1]

When I google, I found [this figure][2] is easy to understand, but the figure description does not make sense for me.

for e.g.
If we take two-fold change;
isn't this the way we calculate LFC;

LFC = Log2 (Foldchange) = Log2(2) =>1

if X axix is in Log2 fold change, shouldn't two-fold change denotes at **± 1**? Why they denote 2 fold change cutoff line in **± 2**?

-----------------
**Below VP; Fold Change greater than 2 times and significant P value less than 0.05**

See this example:

![enter image description here][3]

and this:

![enter image description here][4]

Could someone help me to understand this terminology ( **"Fold change", "Log Fold Change" *(LFC)*, and "Percentage of Change"**. ) easily with an example?

[1]: http://www.nathalievialaneix.eu/doc/pdf/tutorial-rnaseq.pdf
[2]: https://documentation.partek.com/display/FLOWDOC/Volcano+Plot
[3]: https://cdn.xlstat.com/uploads/assets/Tutorials/DifferentialExpression/differentialExpression06.png
[4]: http://chenyuan.date/images/Volcano-plot-unnamed-chunk-4-1.png

Answer: All the figures you posted above shows same thing which is Log2(FoldChange). Although 1st figure x axis labeled "Fold Change" it is Log2(FoldChange). Without log you cannot have negative values (It means down regulated genes).

One of the reasons to use log2(FoldChange) and not FoldChange only is to use same magnitude for UP and DOWN regulated genes. See the example below.

Let's say you have 2 genes one **GeneA** and **GeneB**.

**GeneA** has FPKM (or RPKM) value is **100 in control** and **value 200 in treatment**.
`FoldChange for GeneA is : (treatment/control = 200 / 100 = 2 )`.

**GeneB** has FPKM (or RPKM) value is **200 in control** and value **100 in treatment**.
`FoldChange for GeneB is : (treatment/control = 100 / 200 = 1/2)`.

As you can see from FoldChange values, **GeneA is 2 fold UP** while **GeneB is 1/2 (0.05) fold down**. Comparison of value 2 and 0.05 are less interpretable to indicate UP and DOWN gene regulation. However, normalising both UP and DOWN on log2 make UP and DOWN comparisons easy. For e.g. Log2FC for GeneA is 1 and GeneB is -1 which makes fold change comparable and so easy to plot figures (volcano , heatmap etc. ).

In most RNAseq people use 1.5 or 2 foldChange cutoff. It means Log2FC 0.6 and 1 respectively.

In short, log2FC value >= 0.6 & =< -0.6 suggests 1.5 fold cutoff and log2FC value >= 1 & =< -1 suggests 2 fold cutoff.
Question: Dear Biostars,

Could I please know, If there is some way to view the reads in BAM file according to user defined Insert Size range. I would like to visualize all the reads falling between certain Insert Size ranges, not all the reads.

Your ideas and resources would help me a lot. Thanks in advance.

~Prakki Rama.

Answer: This is a hack but it might be good enough? Parse the input bam to extract alignments within a given *template* size range, here between 0 and 200:

samtools view -h read.pe.bam \
| awk '{if (($0 ~ /^@/) || (sqrt($9^2) > 0 && sqrt($9^2) < 200)) {print $0}}' \
| samtools view -Sb - > insert.bam &&
samtools index insert.bam

Now load insert.bam in your favourite browser.
Question: Really simple question but for some reason I cannot find a simple answer to this.

We are moving from SGE to SLURM on our HPC. On SGE, we often submit a job with a range of threads, then inside the job, read the number of threads we got and use that to configure multithreaded programs. Example:

$ cat test.sh
#!/bin/bash

echo "got $NSLOTS cpu slots"

$ qsub -pe threaded 4-16 test.sh
Your job 5597119 ("test.sh") has been submitted

$ cat test.sh.o5597119
got 15 cpu slots

What is the SLURM equivalent of this?

I have tried a million variations of `sbatch test.sh --cpus-per-task x --ntasks y` across different queues, but nothing seems to work. Its also not clear how to request a range of threads, nor is it clear how to figure out how many you got; I've seen examples online that use `$SLURM_NTASKS` and `$SLURM_CPUS_PER_TASK` but neither of those variables exist inside any of my jobs. The SLURM [docs](https://slurm.schedmd.com/mc_support.html) have been pretty [unhelpful](https://slurm.schedmd.com/sbatch.html).

Answer: You can't use a range of threads, you must state how many you need. The program you run in the script only knows how many threads you've allocated if you tell it. It is completely possible to lie to slurm and say:

#SBATCH -c 1
bowtie2 -p 20 ...

Jobs won't be killed for using too much CPU (it's often impossible to predict exactly how many will get used), at least on our cluster, rather what you specify with `SBATCH -c X` is removed from the list of allocatable resources.

Make sure to specify `#SBATCH --ntasks-per-node=1` to prevent tasks from getting split across nodes (unless your processes can actually handle this...few can in bioinformatics).
Question: I'm having problems with **quality score recalibration** using **GATK** as most tutorials/examples use an old version of GATK which has different syntax and arguments to the current version.

I have run this to get a recalibration table, which I think I need to get the recalibrated BAM files:

gatk-4.1.2.0/gatk BaseRecalibrator \
-reference human_genome/hg38_masked.fa \
--known-sites human_genome/All_20180418.vcf.gz \
--input S01_BAM_file.bam \
--output recal.table

This works fine, and I could run this for all of my samples.

But when I then do the next step to get the **recalibrated BAM files** I can't work out what to do.

*The old way of doing it was:

java –jar GenomeAnalysisTK.jar –T PrintReads \
–R human.fasta \
–I realigned.bam \
–BQSR recal.table \
–o recal.bam*

I wrote:

gatk-4.1.2.0/gatk PrintReads \
-R data_R3/human_genome/hg38_masked.fa \
--input S01_BAM_file.bam \
-**BQSR** recal.table \
--output quality_score_recalibrated_S01.bam

But this makes an error that says:
**"A USER ERROR has occurred: B is not a recognized option".**
And I can't find an alternative argument for -BQSR when I look at the PrintReads help manual.

**Does anyone know the current way to run this quality score recalibration to get recalibrated BAM files?**

Answer: you should use ApplyBQSR :

https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_bqsr_ApplyBQSR.php

gatk BaseRecalibrator \
-I input.bam \
-R reference.fasta \
--known-sites sites_of_variation.vcf \
--known-sites another/optional/setOfSitesToMask.vcf \
-O recal_data.table

then

gatk ApplyBQSR \
-R reference.fasta \
-I input.bam \
--bqsr-recal-file recal_data.table \
-O output.bam
Question: DESeq2 gives logFC values based on fitted model. We can also get the normalized counts for each sample based on the raw counts and normalization factors, thus we can calculate the logFC ourselves using normalized counts. The question is how to report the final result. I prefer the logFC+pval/FDR from the fitted model while some colleagues prefer the self-calculated logFC and model-based pval/FDR (usually these logFC values are more dramatic than that from the fitted model). I feel it is a bit weird to combine different sources of results. What do you usually do? Will you combine FC values from one method and p values from another? Just curious to know which is a common practice. Thanks.

Answer: Please ignore your colleagues, they're doing it wrong.

One of the benefits of DESeq2 is that it shrinks fold-changes (toward a Gaussian prior if I remember correctly), which should increase their reliability. This is discussed and tested in [the preprint][1], which I advise everyone that uses DESeq2 to read (for a paper on a statistics package, it's surprisingly accessible).

So just report the shrunken fold-changes and adjusted p-values reported by DESeq2 and politely hand your colleagues a copy of the preprint.

[1]: http://biorxiv.org/content/early/2014/02/19/002832
Question: I always assumed the designations of "forward" and "reverse" strand on DNA were chosen arbitrarily (i.e. the strands could just as easily have been named the other way around).

However, a quick browse of the chromosome maps for Human on the NCBI map viewer seems to imply that in all cases, base numbering starts on the side closest to the centromere. Assuming that this numbering applies to the "forward" strands of the reference genome, does this mean that, for humans at least, there is a policy whereby the forward strand is always the one which contains the centromere at the lowest numerical index?

Also, if this is the case, does anybody know if this convention would apply outside humans?

Thanks for your time.

Answer: > ...does this mean that, for humans at least, there is a policy whereby the forward strand is always the one which contains the centromere at the lowest numerical index?

Yes, in humans the centromere is always closest to the start coordinate of each chromosome, since human chromosomes are oriented based on their p (for "petit") and q arms. The p and q designation comes from the size of the chromosomes arms observed classically during mitotic cytogenetic preparations.

> Also, if this is the case, does anybody know if this convention would apply outside humans?

No, this convention does not apply generally outside of humans. For example, the yeast Saccharomyces cerevisae has sixteen chromosomes, 9 of which have the centromere closer to the start of the chromosome than the end.

```
$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D sacCer1 -e 'select sgdOther.chrom as "chrom.", size as "chrom. size", sgdOther.chromStart as "left arm", size-chromStart as "right arm", chromStart/(size-chromStart) as "left/right ratio" from sgdOther, chromInfo where name like "CEN%" and sgdOther.chrom=chromInfo.chrom order by size/(size-chromStart);'
+--------+-------------+----------+-----------+------------------+
| chrom. | chrom. size | left arm | right arm | left/right ratio |
+--------+-------------+----------+-----------+------------------+
| chr12 | 1078173 | 150826 | 927347 | 0.1626 |
| chr8 | 562639 | 105578 | 457061 | 0.2310 |
| chr5 | 576869 | 151959 | 424910 | 0.3576 |
| chr13 | 924430 | 268030 | 656400 | 0.4083 |
| chr2 | 813136 | 238167 | 574969 | 0.4142 |
| chr4 | 1531914 | 449706 | 1082208 | 0.4155 |
| chr15 | 1091285 | 326583 | 764702 | 0.4271 |
| chr3 | 316613 | 114378 | 202235 | 0.5656 |
| chr7 | 1090944 | 496920 | 594024 | 0.8365 |
| chr6 | 270148 | 148503 | 121645 | 1.2208 |
| chr10 | 745446 | 436001 | 309445 | 1.4090 |
| chr16 | 948060 | 555951 | 392109 | 1.4178 |
| chr1 | 230208 | 151456 | 78752 | 1.9232 |
| chr11 | 666445 | 439773 | 226672 | 1.9401 |
| chr14 | 784328 | 628756 | 155572 | 4.0416 |
| chr9 | 439885 | 355625 | 84260 | 4.2206 |
+--------+-------------+----------+-----------+------------------+
```
Question: **Background**

I have paired-end RNA-seq reads from a drug-treatment experiment, with < 15 million mapped reads in many samples (too few reads) and large variability in mapped reads across biological replicates. Differential expression and splicing analysis on these samples indicate that statistical power in my tests could be improved if I had better sequencing depth, and I have remaining RNA from these samples available to re-sequence.

**The questions**

Is it analytically and statistically tractable to re-sequence the same samples and control for potential artifacts in the combined data?

What would be the best workflow for merging data from these two RNA-seq runs? I would guess that it's best to keep the runs separate until the counts have been summarized. Then I can carry out PCA to visually inspect the gross extent of artifact in the different runs before merging the counts for statistical analyses.

Beyond gross visual inspection of PC's, what sorts of quality control steps could I take if I identify a strong batch effect between the different sequencing runs? Would software like svaseq or combat be appropriate here if I do identify a batch effect? If so, would it be best to remove the batch effect in the samples before combining the count data?

Answer: As you suggested (and confirmed by genomax2) it's probably the best to check using PCA if your two runs result in approximately the same result.

But as soon as you have determined it's okay I would suggest to merge your bam files, and repeat the counting before you do your final analysis. That would minimize your chance of errors.

Furthermore, in case you are using a two-step alignment (e.g. using STAR) it might be advantageous to merge the fastq files across runs and repeat the alignment.
Question: I don't know how stupid of a question this is so please bear with me. I have a certain number of EST cDNA sequences that I used for miRNA prediction. The results are cDNA sequences so for the secondary structure prediction of precursor miRNA, should I use the RNA folding form or the DNA folding form on Mfold because it gives both options and I am confused on which should I use. DNA folding because my sequences are cDNA sequences or RNA folding form because I am trying to predict pre-miRNA sequences after all. Please advice.

Answer: I think you should use RNA parameters because miRNAs are indeed RNA, that doesn't depend on how the sequence was acquired.
Question: Hi all.

I think that p-value is one of the most greatest way of measuring degree of observed data.

However, BLAST doesn't use p-value but E-value.

Why the BLAST use e-value for interpreting sequence data instead of p-value?

Is there any logical reason to use E-value for BLAST? If so, could you tell me the detail reason?

Answer: Quote from the BLAST help (http://www.ncbi.nlm.nih.gov/BLAST/tutorial/Altschul-1.html#head4 ):

"The BLAST programs report E-value rather than P-values because it is easier to understand the difference between, for example, E-value of 5 and 10 than P-values of 0.993 and 0.99995. However, when E < 0.01, P-values and E-value are nearly identical."

Important to note that P value of the BLAST is not the same thing than a P-value of a t-test.
Question: Hi,

I was wondering if there is a way I can generate junctions.bed file from hisat2 just like the one in tophat.

Hisat2 documentation says use extract_splice_site.py python script but that gets splice site from gtf file, is there a way to know splice sites and no of reads associated with those sites using hisat2?

Thanks

Answer: Easiest solution will be to use [featureCounts][1] and counts reads from exon-exon junctions by providing your aligned BAM files. Check manual [page 35][2] for detailed information.

**-J (juncCounts)** Count the number of reads supporting each exon-exon junction.
Junctions are identified from those exon-spanning reads
(containing ‘N’ in CIGAR string) in input data. The output
result includes names of primary and secondary genes that
overlap at least one of the two splice sites of a junction. Only
one primary gene is reported, but there might be more than
one secondary gene reported. Secondary genes do not overlap
more splice sites than the primary gene. When the primary
and secondary genes overlap same number of splice sites, the
gene with the smallest leftmost base position is selected as
the primary gene. Also included in the output result are the
position information for the left splice site (‘Site1’) and the
right splice site (‘Site2’) of a junction. These include chromosome
name, coordinate and strand of the splice site. In
the last columns of the output, number of supporting reads is
provided for each junction for each library

[1]: http://bioinf.wehi.edu.au/featureCounts/
[2]: http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf
Question: Hello everyone,

I have a bunch of TARGET samples but the associated annotation lacks the tissue information i.e. whether it is a Primary tumor or Recurrent etc. I know the value for a couple codes like `01` stands for Primary solid tumors, but I want to define the tissue type for every code for future references. Does anybody know where I can find the TARGET code table to interpret the codes?

TARGET_10 (Blood Derived Normal)
TARGET_14 (Bone Marrow Normal)
TARGET_11 (Solid Tissue Normal)
TARGET_01 (Primary Solid Tumor)
TARGET_09 (Primary Blood Derived Cancer - Bone Marrow)
TARGET_03 (Primary Blood Derived Cancer - Peripheral Blood)
TARGET_04 (Recurrent Blood Derived Cancer - Bone Marrow)
TARGET_02 (Recurrent Solid Tumor)
TARGET_40 (Recurrent Blood Derived Cancer - Peripheral Blood)
TARGET_06 (Metastatic)
TARGET_41 (Blood Derived Cancer - Bone Marrow, Post-treatment)
TARGET_42 (Blood Derived Cancer - Peripheral Blood, Post-treatment)
TARGET_13 (EBV Immortalized Normal)
TARGET_15 (sample type 15)
TARGET_60 (Primary Xenograft Tissue)
TARGET_50 (Cell Lines)
TARGET_61 (Cell Line Derived Xenograft Tissue)
TARGET_20 (Control Analyte)

Thanks!

**UPDATE**: I have updated the code table based on @khericlim's response.

Answer: Hey Komal, how are you? Say hi to Pichai!

The TCGA one was gone, but this should help.

https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes
Question: In Jellyfish's original paper(Marcais et al, Bioinformatics 2011), it says that "counting kmers up to 31 bases in length".

But on the other hand, in other papers such as KMC2 (Deorowicz et al, Bioinformatics 2015), Jellyfish2 is used to count kmer size = 55.

So does Jellyfish support k > 31 now? and How it is supported?(i.e. why in original version it is not supported but now it is supported)?

Thanks.

Answer: Probably because they originally used 64bits per kmer, which can store up to 32 bases in 2bit (less one for whatever reason, probably to signal reverse complementation or variable-lengths).

With an update, this man-made constraint was removed.
Question: hello

I have a question while i am reading a paper that is related to a gene fusion discovery.

below is some of the paper's content.

> A paired end read will be a member of multiple valid clusters as a result of homology between genes, exon redundancy between transcripts, and the multiplicity of valid clusterings. Fusion events are rare when compared with the event of detecting a discordant paired end read given the existence of a fusion event that would generate that read. **Thus we seek an assignment of each paired end read to a single
fusion event (valid cluster) that minimizes the number of fusion events**.

(Valid cluster represents a potential fusion event implied by a set of discordant alignments)

i don't understand the bold text. I think the verb "minimizes" should be changed to "maximizes " because

It maybe increase the number of gene fusion by calculating the single paired end read to a single gene fusion.

I hope you understand my basic questions. I'm beginner in this field.

Anyway I am looking forward to your reply.

Thank you.

Answer: The problem that the authors are trying to address is the problem of false positive results in the setting of reads that can map to multiple locations. They are suggesting to use a heuristic that takes all the data available and seeks a minimal set of fusions from the set of all potential fusions. So, I believe "minimizes" is what the authors meant.
Question: Hello. I have various .bed and .tped files (+ corresponding .fam/.tfam files) that I am exporting in vcf format from plink as follows:

./plink --tfile examplesnps --chr-set 39 --out exported_vcf --recode vcf

However, I don't know how to get an output that has the sample IDs attached to the genotype data. Does anyone know how to include the actual sample names (i.e., the individual IDs) in the outputs from plink?

Thanks in advance!

Answer: Have you tried `--recode vcf-iid`?
Question: I am a completely rookie on Bioinformatics, so please bear with me and use simple language (I am a computer scientist) :-)

How can we use k-mers to find out if a gene is similar to our query string?

For example: We have a reference gene r = `ACAAGTC`, and a query string q = `CATGT`. For sequence alignment we could get the two possible solutions: (see photo on link)

![enter image description here][1]

But if we use k-mer with k=2 (reason for 2 is that I tried using k=3 and got 1 match) we could split q into the set = {CA, AT, TG, GT} we can get 2 matches.

I am confused on how to use k-mer for queries. I am guessing that it would also need to know where in the reference string it is found, since order matter. But most importantly, **Why** can we use K-mers? That's maybe my biggest question mark.

[1]: https://images2.imgbox.com/bc/62/X0aa6MzH_o.png

Answer: You wouldn't generally use kmers this short as the odds of finding spurious matches becomes massive, but I accept that this is a toy example.

To my knowledge, we don't really approach a sequence alignment with kmers, though you might use them to seed the alignment (in a **local** alignment, small subregions of alignment are found and then slowly expanded).

You can however compute useful statistics using kmers about how similar 2 sequences are when they are inefficient to align (e.g. a whole genome). It's not uncommon to look at sequence 'distances' in this way. You might find this code interesting: https://github.com/jrjhealey/bioinfo-tools/blob/master/StringComparisons.py#L158-L181

Those are a few of the most common kmer string/sequence comparison methods implemented in python.

You can use kmers to get a rough approximation of 2 large sequences in a very computationally efficient manner, where size or number of sequences becomes prohibitive, with the idea being that the more similar 2 sequences are, the more likely you are to randomly select the same kmer from each of them and the more frequently they will probably occur.

Not sure if that answers your question or not?
Question: Dear all,

I have FASTQ files and on start of my read I have 7 nucleotides tag - I would like to extract reads with this specific tag.

I would like to search in first 15 nucleotides of my reads, if match - extract this read to new fastq files.

Thank you for any ideas or help.

Answer: You can use Heng Li's bioawk:

To check if the tag is part of the first 15 bases

bioawk -c fastx 'substr($seq,0,15) ~ /$TAG/ { print }' reads.fq.gz

To match the first 7 bases to your tag,

bioawk -c fastx 'substr($seq,0,7) == $TAG { print }' reads.fq.gz