cannin/biostars_qa huggingface dataset:
Question: Dear all,

I have paired-end fastq data generated with Illumina bcl2fastqv2.19 & sequenced on a Novaseq.The i5index is 7bp long, the i7 8bp long

R1.fastq.gz contains R1 101bp reads:

@A00154:125:HGKTMDMXX:1:1101:10420:1000 1:N:0:AACTGAGG+ATGCGTC
CTGGCCGTCTCAGCCGAGAAGCCGAGGATTGAATGGGCATGGAGACTGAACTACCCCTCTCACCTTTAGAGGTGGCTCCTCCAAGTCGGGGTTGACGCCCG
+
FFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

R2.fastq.gz contains 6bp UMI sequence

@A00154:125:HGKTMDMXX:1:1101:10420:1000 2:N:0:AACTGAGG+ATGCGTC
GCGCGT
+
FFFFFF

R3.fastq.gz contains R2 101bp reads:

@A00154:125:HGKTMDMXX:1:1101:10420:1000 3:N:0:AACTGAGG+ATGCGTC
CTTCATAGGCCACAAAAAGCCCATATATCAGTGTCATCCACTAAGCCTCAGACACTGCAGCACGGGCAGCGGCAGTGCCAGCTTCGCCCACACTGCCCCTC
+
FFFFFFFFFFFFFFFFFFFFFF:FF:FFF:FFFFFF:FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

In a downstream analysis I want to use UMI-tools for deduplication. However for that I need the UMI be part of the read name.
@Instrument:RunID:FlowCellID:Lane:Tile:X:Y:**UMI** ReadNum:FilterFlag:0:IndexSequence or SampleNumber

There are tools to add a UMI to the read name when the UMI is present in the read itself. But in my case, the UMI is in a seperate fastq. How could this be achieved?

Answer: An `awk` solution:

$ awk -v FS="\t" -v OFS="\t" 'NR==FNR {split($1, id, " "); umi[id[1]]=$2; next;} {split($1, id, " "); $1=id[1]":"umi[id[1]]" "id[2]; print $0}' <(zcat R2.fastq.gz|paste - - - -) <(zcat R1.fastq.gz|paste - - - -)|tr "\t" "\n"|bgzip -c > R1_umi.fastq.gz

$ awk -v FS="\t" -v OFS="\t" 'NR==FNR {split($1, id, " "); umi[id[1]]=$2; next;} {split($1, id, " "); $1=id[1]":"umi[id[1]]" "id[2]; print $0}' <(zcat R2.fastq.gz|paste - - - -) <(zcat R3.fastq.gz|paste - - - -)|tr "\t" "\n"|bgzip -c > R3_umi.fastq.gz

The `fastq.gz` files get uncompressed by `zcat` and the 4 line belonging to a read get tab delimited by `paste`.

`awk` saves the id and the umi in a list, where the key is the header until the first white space, and the value is the umi code.

If the second `fastq` file is read, we append the umi to the id and print out the line. Then the tabs are reverted to new lines by `tr` and the file get compressed using `bgzip`.

fin swimmer
Question: Hello everyone,

I know that affymetrix probeset ids are portable across platforms: the same probeset id `1090_f_at` always refers the same gene (probe set). Is it the same for probe ids?

Executing, as discussed in the "affy" package vignette, the code:

> pm(Dilution, "1090_f_at")

on the chip HG_U95Av2, for the probeset 1090_f_at are returned its pm probes 1, 2, .. 16

```
20A 20B 10A 10B
1090_f_at1 115.0 74.0 94.0 61.0
1090_f_at2 129.3 80.3 108.0 70.3
..
1090_f_at16 262.0 157.0 194.8 131.5
```

Executing the same code on a different chip containing, of course, the same probeset `1090_f_at`, are its pm probes returned in the same order, so that their ids

```
1090_f_at1,
1090_f_at2,
..
1090_f_at16
```

are consistent across chips?

I suppose so. Could anyone confirm it? How could I check it?

Thank you

Answer: It may depend on the source of the file you are using to obtain probe ids. Probe IDs in files from Affymetrix are integers and are associated with a specific feature (x, y position) on the microarray. For clarity this is different than a probeset id. A probeset id refers to a collection of probes (you seem to know this difference but I see people use the terms probe and probeset interchangeably in casual conversation which leads to confusion...)

Probe IDs for Affymetrix arrays (in files released from Affymetrix, this may differ depending on the creator of your analysis package) are often calculated based on a formula and convertible from the probe id to an x y coordinate on the microarray.

Therefore the only reason a probe_id should be the same between two different arrays is if the probe is the same height and width, and the probe is at the same x, y coordinate on the microarray. (This is rare).

You can download lots of information from the Affymetrix website. If you have any further questions just ask. I'm pretty familiar with it.
Question: Hi,

I'm working on bacterial sulfur metabolism and I would like to be able to extract as much information as I can from the different pathways linked to this metabolism (at least all the biochemical reactions and the enzymes that catalyze these reactions and if possible, the genes coding for these enzymes). Here is what I've done so far:

I've obtained a free license to download the [MetaCyc][1] database in Biopax level 2 and 3 format so I try to extract the information from these files.

I've used the R package [rBiopaxParse][2]r to read and import the Biopax in R but now **I'm quite stuck because I don't really know how I can proceed to extract a pathway** ( by pathway I mean the pathway level but also all the lower levels of information like nested pathways, reactions, molecules and proteins interactions,...).

My aim would be to get these pathways as individual Biopax files and the covert them in sif format thanks to the R package [paxtoolsr][3] and then visualise my pathways as a network and manipulate it with graphiz or [RCytoscape][4] packages.

I don't know if it's the best strategy to be able to visualise and manipulate the whole metabolic pathways but since I don't have any knowledge in XML like langages I try to use R to get rid of this format^^. What could I do to be able to extract the needed information to get my metabolic pathways networks (all the biochemical reactions and the enzymes that catalyze these reactions)?

Best regards,
[1]: http://metacyc.org/
[2]: https://bioconductor.org/packages/release/bioc/html/rBiopaxParser.html
[3]: http://bioconductor.org/packages/release/bioc/html/paxtoolsr.html
[4]: https://www.bioconductor.org/packages/release/bioc/html/RCytoscape.html

Answer: Try this just using paxtoolsr with the BioPAX Level 3 file you have:

# Update paxtoolsr

First, update to the paxtoolsr development version (I just updated a few things):

setRepositories(ind=1:6)
options(repos="http://cran.rstudio.com/")
if(!require(devtools)) { install.packages("devtools") }
library(devtools)
install_github("BioPAX/paxtoolsr")

# Extract Pathways

library(paxtoolsr)

# An example with the sample BioPAX file in the paxtoolsr package
exampleFileInPaxtoolsr <- system.file("extdata", "REACT_12034-3.owl", package="paxtoolsr")
sifnx <- toSifnx(exampleFileInPaxtoolsr, "output.txt", "uniprot")

# Not all rows represented, but that's because not every row has a pathway listed
rowIndiciesForPathways <- splitSifnxByPathway(sifnx$edges)

# A pathway extracted
bmp <- sifnx$edges[rowIndiciesForPathways$`Signaling by BMP`, ]

# If you prefer a data.frame over a data.table (data.table is used for file reading speed), then do this:
library(data.table)
class(bmp) # Should be "data.table" "data.frame"
setDF(bmp)
class(bmp) # Should be "data.frame"

# Plot in R

# For simple plotting in R
g <- loadSifInIgraph(bmp)
plot(g)
Question: **Edit 1:** The following is Y/Auto depth versus age in a cohort. Individuals in this cohort are not known to be cancer patients.

![enter image description here][1]

**Original question:**

In WGS/WES studies I read and studies I participated in analyzing, I notice that the average sequencing depth on Y chromosome vary a lot more than X chromosome, and the depth is mostly less than the average depth on autosomes (see the figure for an example; from [Turro et al. 2020][2]). Most individuals that got sampled do not have diseases involving large deletion of Y chromosome.

What causes the variability in the sequencing depth on Y chromosome? Why is the average depth on Y almost always less than that on autosomes?

![enter image description here][3]

[1]: /media/images/50b4a58e-0c41-4919-890c-49cd71bc
[2]: https://www.nature.com/articles/s41586-020-2434-2.pdf
[3]: /media/images/835d5568-c6e1-4f3f-a0b0-2836db5b

Answer: Hi Samuel,

I would guess it is the age related "Loss of the Y chromosome" you see in the data, which is a known effect.
See:
[Loss of Chromosome Y and Its Potential Applications as Biomarker in Health and Forensic Sciences][1]
for example.

"Loss of chromosome Y (LOY) is a mosaic aneuploidy that can be detected mainly in blood samples of male individuals. Usually, LOY occurrence increases with chronological age in healthy men."

"The mechanism by which LOY occurs is yet to be elucidated. It is not known whether the Y chromosome is lost as a whole or if it is a gradual deletion process during the lifespan."

We see that quite a lot in our WES data of older men.

[1]: https://www.karger.com/Article/FullText/508564
Question: Dear everyone,

I have been asked to comment on an experiment design that involves WES or WGS of cancer cell lines which lack matched normal.

I know that this design is far from ideal, but I was wondering if there are people who have already stream-lined it (to the extent possible of course). I can think of a couple of common variation filters to keep out the germline variants, but I am not sure how to go beyond.

I am sorry if this is a duplicate question; I couldn't find a related hit.

Thank you,

Noushin

Answer: The bottom line is that without a matched normal, you're just not going to be able to call the somatic status for the vast majority of sites. That said, you can winnow down a list to those you *suspect* are somatic. Some ideas:

Weed out sites with high frequency in the population

If your tumor is very impure, you can take advantage of the fact that the frequencies of somatic variants will be shifted away from 50%/100%
Question: I was checking some of the variants I have found in cancer samples that I have been analysing. I am not much aware of how COSMIC data should be interpreted. So if something is already found in COSMIC, can we consider that variant novel or something significantly important for publication?

Answer: COSMIC does provide a large catalogue of somatic variants from many sequencing studies, although their data can be somewhat messy (they include variants sometimes which are not marked as somatic). Examining COSMIC can be reasonably good evidence at determining whether somatic variants have been seen before. That said, the majority of somatic variants in COSMIC are passenger mutations, and even if a hypothetically true driver mutation has been seen in COSMIC a few times does not mean people have previously implicated it as a driver (you might have to backtrack the study to examine their conclusions though). In fact, rare driver genes (especially <5-10% of samples) are generally poorly characterized in most tumor types (http://www.nature.com/nature/journal/v505/n7484/abs/nature12912.html). COSMIC updates regularly, but the current version may not have included all studies recently performed. So I would suggest that a literature search for very recent studies related to your particular cancer samples would still be warranted. Lastly and perhaps most important, whether your variants are "important" greatly matters on whether they are driver or passenger mutations. A passenger mutation that has not been seen before is not really important. Whether mutations may be drivers can be supported statistically/computationally and/or by experiments. Some support for driver status may be whether your cancer samples support the gene as significantly mutated (e.g. MutsigCV, MUSIC), high in silico predicted score for mutations (e.g. [CHASM, VEST](http://www.cravat.us/)), and for protein coding genes where the mutation occurs in 3D protein structure ([MUPIT](http://mupit.icm.jhu.edu/) maps cancer mutations to 3D protein structure with functional annotations).
Question: I am 34 now. I have no time to prepare for any cet for PhD. I need a PhD for my job and my future.

Please can you help me how to do PhD in bioinformatics, and where to go for distance education in bioinformatics?

Answer: > I am 34 now .i have no time to prepare for any cet for phd

Its never too late. See [this][1]

> ..i need a phd for my job and my future .

No, that cant be true. People have jobs and future even without a PhD.

> please can help me how to do phd in bioinformatics ,and where to go
> for distance education in bioinf

You are missing out the whole idea of a PhD. It's not like any other degree where you opt for some courses, you sit in some exams, you are graded and then you have a certificate. A "PhD" is a process; a journey where you learn by experience and failures. You learn to design your experiments; frame your own questions; develop a sense of problem solving. You fail at times and learn by mistakes and move on. That entire process ; that journey, polish you as a scientist. You develop abilities that you apply in other areas when you done with it. You miss all of this is a part time or distance PhD.

[1]: https://chrisblattman.com/2013/06/12/when-are-you-too-old-for-a-phd/
Question: The invocation of [VarDictJava][1] in paired variant calling mode is

VarDict -b "/path/to/tumor.bam|/path/to/normal.bam"

Now obviously in my CWL I would like the tumour and normal BAM to be separate inputs. However, I'm not sure how to combine these into a single command line argument in the actual CommandLineBinding

[1]: https://github.com/AstraZeneca-NGS/VarDictJava#paired-variant-calling

Answer: Hello @ttmigueltt,

Here's an alternative solution if you will always have both samples; no `InlineJavascriptRequirement` needed:

cwlVersion: v1.0
class: CommandLineTool

inputs:
normal:
type: File
format: edam:format_2572
tumor
type: File
format: edam:format_2572

baseCommand: echo

arguments:
- $(inputs.bam_input.path)|$(inputs.bam_input.path)

outputs: []

$namespaces: { edam: http://edamontology.org/ }
$schemas: [ http://edamontology.org/EDAM_1.20.owl ]
Question: When trying to investigate a 16S rRNA dataset, I often identify several dozen/hundred species/families which are found in higher/lower abundances. I then start doing literature searches to see what they could be doing, where they have been observed before etc.

To me this sounds:

Really selective, only sampling a few papers for each species.

Limiting, as there is no way to do this fully for tens of species.

Incredibly time consuming.

What I'm really looking for is a system which I can put a taxa list into and it'll say "Those ones are all anoxic" or "Those 5 have shown denitrifying ability". I don't know if this could be done with literature mining or where to start this, or if there is a database around in the world which curates data like this...

Any suggestions are appreciated.

Answer: On following Neil's advice, and reading some papers I reverted back to IMG to see what meta-data they collect with their published genomes, and actually, it's quite substantive.

On clicking into the 'Genome Browser' section you get presented with a list of species and some basic data on the sequencing project. But if you navigate to the bottom of the page, there is a selection as such:

This allows you to generate your table with the data you're interested in. A selection as above results in the following output, which can then be exported to excel or tsv format and used however you would like. (Full size available here:http://i.imgur.com/KjMXRDu.png)
Question: Dear lazyweb, here is the new menu for `samtools index`:

```
Usage: samtools index [-bc] [-m INT] [out.index]
Options:
-b Generate BAI-format index for BAM files [default]
-c Generate CSI-format index for BAM files
-m INT Set minimum interval size for CSI indices to 2^INT [14]
```

What is the **CSI** index compared to **BAI**? Should I better use it instead of the 'old' BAI index? What's the benefit?

Answer: It depends on how long your contigs/chromosomes are. The biggest benefit to CSI is that it supports indexing BAM files with contigs longer than 2^29-1 bases. At least in plants this isn't an uncommon requirement.
Question: As the title says, I'm working on variant calling for somatic variant discovery where I have tumour samples, but no normal samples to compare with.

Previously I've been using GATK's UnifiedGenotyper for variant calling, but as this tool is deprecated, I want to switch to the newer GATK variant calling tools. However, I can't decide whether it would be more appropriate to use HaplotypeCaller or Mutect2.

On the one hand, HaplotypeCaller can be run without requiring a normal sample for comparison, but the documentation for HaplotypeCaller states that:
> the algorithms used to calculate variant likelihoods is not well suited to extreme allele frequencies (relative to ploidy) so its use is not recommended for somatic (cancer) variant discovery. For that purpose, use MuTect2 instead.

However, Mutect2 is still in Beta, and although it can be run on tumour samples only, and in the Mutect2 documentation it says:
> Tumor-only variant calling is possible but it is NOT supported and we will not answer any questions about it until it becomes a supported feature.

So neither tool is ideal for my purposes, but I'm having trouble deciding which will be the most applicable to my data.
Any suggestions or advice would be greatly appreciated.

Answer: Some resources that address your question more broadly:

- Evaluating Variant Calling Tools for Non-Matched Next-Generation Sequencing Data: https://www.nature.com/articles/srep43169
- https://www.biostars.org/p/207536/
Question: Hi all,

Can someone help me understand the RSeQC Output from infer_experiment.py?

So this is the output:

This is PairEnd Data
Fraction of reads failed to determine: 0.0560
Fraction of reads explained by "1++,1--,2+-,2-+": 0.0192
Fraction of reads explained by "1+-,1-+,2++,2--": 0.9247

So it's stranded but is it fr-firststrand or fr-secondstrand?
I do not understand the help given [here][1]:

> For pair-end RNA-seq, there are two different ways to strand reads
> (such as Illumina ScriptSeq protocol):
>
> 1++,1–,2+-,2-+
>
> read1 mapped to ‘+’ strand indicates parental gene on ‘+’ strand
>
> read1 mapped to ‘-‘ strand indicates parental gene on ‘-‘ strand
>
> read2 mapped to ‘+’ strand indicates parental gene on ‘-‘ strand
>
> read2 mapped to ‘-‘ strand indicates parental gene on ‘+’ strand
>
> 1+-,1-+,2++,2–
>
> read1 mapped to ‘+’ strand indicates parental gene on ‘-‘ strand
>
> read1 mapped to ‘-‘ strand indicates parental gene on ‘+’ strand
>
> read2 mapped to ‘+’ strand indicates parental gene on ‘+’ strand
>
> read2 mapped to ‘-‘ strand indicates parental gene on ‘-‘ strand

Thanks for your help!

[1]: http://rseqc.sourceforge.net

Answer: It means you have a standard (dUTP-based) strand-specific library. If you want to use featureCounts, you'll want the `-s 2` setting. For HTSeq-count it's `--strand reverse`.
Question: I have bam files that I have split into unmapped, uniquely mapped and multimapped bam files from HISAT2 alignment. I am trying to merge the sorted unmapped and uniquely mapped bam files with the command

samtools merge A1_merged.bam -b A1_unmapped.bam A1_unique.bam

and am getting an error

[E::hts_open_format] Failed to open file BAMT
samtools merge: fail to open "BAMT": No such file or directory

the square actually has a question mark in it but that didn't copy across. Anyone have any ideas on why this error is coming up and how to fix it?

Answer: -b FILE List of input BAM filenames, one per line [null]

`-b` would be necessary if the file names were in a separate text file.

It is simply

samtools merge A1_merged.bam A1_unmapped.bam A1_unique.bam
Question: Hello, this is a very basic question but I was wondering if someone could help me understand if I've used the correct GTF file and Fasta file for the mouse genome indexing. I got the relevant Fasta file and GTF file from ensembl:
GTF:ftp.ensembl.org/pub/release-103/gtf/mus_musculus/Mus_musculus.GRCm39.103.gtf.gz
Fasta:ftp.ensembl.org/pub/release103/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz

Or shell I use Mus_musculus.GRCm39.dna.toplevel.fa.gz for fasta to make Generating genome indexes in STAR?
STAR --runMode genomeGenerate --runThreadN 8 --genomeDir index_reference --genomeFastaFiles Mus_musculus.GRCm39.dna.primary_assembly.fa --sjdbGTFfile Mus_musculus.GRCm39.103.gtf

Thank you for your help!

Answer: The genome construction step line looks good. Also consider including the `--sjdbOverhang` parameter. While the default value of 99 is usually fine, as stated in STAR's manual this should be chosen according to your maximum read length in the dataset by subtracting 1. So, if you have 100 nt reads, `--sjdbOverhang 99` is fine.

If you performed trimming and has variable read lengths, then you choose this value according to `max(readLength-1)`. So if you have reads varying from 50-120, you'd have `--sjdbOverhang 119` as optimal parameter.
Question: Which version control system do you prefer to use in your projects and which one is more commonly used. I tried using SVN and CVS but never got the hang of them. I currently use git for most of my projects. Which one is more common among bioinformaticians?

Answer: Having made the progression from CVS to SVN to git, I have to say I prefer the latter, although it's taken me a while to wrap my head around it where I always felt quite comfortable with CVS/SVN. Using github has helped a lot, but as I found out yesterday the web based merging tools leave a little to be desired.

I follow dozens of coders and bioinformatics projects on github, you can't deny that the social element makes the system far more valuable than just having an SVN repository tucked away somewhere. I spend a good hour yesterday poring over other peoples code there, it's a great resource, and certainly for NGS tools seems to be a popular repository, backed by a popular SCM.
Question: I work with RNA-seq data and have found a few deferentially expressed genes across particular tissue sample. Now I have been instructed to work with GTEX data to see the deferentially expressed genes across different tissue samples.

Now to go green with GTEX data set, I first don't understand their sample codes like

GTEX-N7MS-0007-SM-2D7W1 for which tissue?

I tried to search the bar codes for GTEX but haven't found any. Please can anyone give some idea for decoding the GTEX bar codes and also performing such analysis. I am sorry if this question is silly because I am a complete naive in the field of NGS.

Answer: In [the download section][1], where it says "A de-identified, open access version of the sample annotations available in dbGaP.", you should find a file called `GTEx_Data_V4_Annotations_SampleAttributesDS.txt`, containing the annotation of each sample. For example GTEX-N7MS-0007-SM-2D7W1 is from Whole Blood.

[1]: http://www.gtexportal.org/home/datasets2