cannin/biostars_qa huggingface dataset:
Question: Hello folks,

I happen to have a small problem which seemed to be trivial at first, but keeps me busy for a while now already. Maybe you can help anew, because most Biostars posts on this issue are quite outdated:

Problem:
-------------------

I need to split 615M paired reads currently in two FastQ files (forward and reversed reads separate) into two file pairs with 308M reads.

Solution attempt A:
-------------------

I unsuccessfully tried to use line count based tools like `split` or `awk`, but the `wc -l` does not even count them correctly (line counts differ), although they are processed nicely by several tools. By googling, I found that a possible reason for this behavior might be that newline characters may occur in the quality scores, which thus hinder the core utilities.

Solution attempt B:
-------------------
I tried to use my `BBTools`, my favorite Swiss army knife of everything sequence related, but...

bbmap/reformat.sh in=... in2=... out=... out2=... reads=308000000
bbmap/reformat.sh in=... in2=... out=... out2=... skipreads=308000000

resulted in

Input is being processed as paired
Input: 615307122 reads 91326404105 bases
Output: 615307122 reads (100.00%) 91326404105 bases (100.00%)

respectively

Input is being processed as paired
Input: 615307122 reads 91326404105 bases
Output: 0 reads (0.00%) 0 bases (0.00%)

for the second command, such that I got a copy of my data in the first case and an empty file in the second.

Solution attempt C:
-------------------

Here on Biostars, somebody had suggested `famas` for a similar problem, so I installed and tried it:

famas --in=... --in2=... --out=...XXXXXX.fq.gz --out2=...XXXXXX.fq.gz -x 308000000

However, the program flooded to output directory with thousands of subfiles files (instead of the actually needed two files each) until the file system couldn't cope with the number of open files anymore and ran out of file descriptors.

ERROR(famas.c|open_output_one:1056): Couldn't open =...compressed.065534.fq.gz
ERROR(famas.c|main:1163): Couldn't open output files. Exiting...

*In case somebody is wondering why I mixed the long and short notation here: the parameter --split-every=308000000 was not recognized (unknown parameter).*

Since it took me a while to clean that mess up again, I am somewhat reluctant to try out more. Does sombeody have any ideas where I screwed up the commands or has suggestions which tools work better?

Thanks a lot for reading, pondering and help!

Matthias

Answer: with [seqkit][1] from the manual:

$ seqkit split2 -1 reads_1.fq.gz -2 reads_2.fq.gz -p 2 -O out

out is output directory, reads_1 is R1 and reads_2 is R2. This function is to break the file into two parts. try `-s` option for sequence number based split.

[1]: https://bioinf.shenwei.me/seqkit/download/
Question: Has anyone had any luck making calls to the UCSC genome API for DNA sequences from the T2T CHM13v2.0 (hs1) assembly? I can easily make calls for hg38 (e.g. api.genome.ucsc.edu/getData/sequence?genome=hg38;chrom=chr1;start=4321;end=5678), however making a call to the hs1 genome doesn't seem to be working. I've checked the list of available UCSC genomes (found here https://api.genome.ucsc.edu/list/ucscGenomes) and hs1 seems to be a valid genome to query, however I am only met with "Bad Request" status messages whenever I try to use it.

I'm sure that I'm just missing something small but some help would be greatly appreciated.

Thanks in advance!

Answer: I couldn't access hs1 too but there is another hub for T2T CHM13v2.0

Listing chromosomes

❯ curl 'api.genome.ucsc.edu/list/chromosomes?hubUrl=https://hgdownload.soe.ucsc.edu/hubs/GCA/009/914/755/GCA_009914755.4/hub.txt;genome=GCA_009914755.4'

Response:

{ "downloadTime": "2022:12:23T22:40:38Z", "downloadTimeStamp": 1671835238, "hubUrl": "https:\/\/hgdownload.soe.ucsc.edu\/hubs\/GCA\/009\/914\/755\/GCA_009914755.4\/hub.txt", "genome": "GCA_009914755.4", "chromCount": 25, "chromosomes": { "CP068254.1": 16569, "CP068257.2": 45090682, "CP068256.2": 51324926, "CP068259.2": 61707364, "CP086569.2": 62460029, "CP068258.2": 66210255, "CP068260.2": 80542538, "CP068261.2": 84276897, "CP068262.2": 96330374, "CP068263.2": 99753195, "CP068264.2": 101161492, "CP068265.2": 113566686, "CP068266.2": 133324548, "CP068268.2": 134758134, "CP068267.2": 135127769, "CP068270.2": 146259331, "CP068269.2": 150617247, "CP068255.2": 154259566, "CP068271.2": 160567428, "CP068272.2": 172126628, "CP068273.2": 182045439, "CP068274.2": 193574945, "CP068275.2": 201105948, "CP068276.2": 242696752, "CP068277.2": 248387328} } %

get sequence:

❯ curl 'api.genome.ucsc.edu/getData/sequence?hubUrl=https://hgdownload.soe.ucsc.edu/hubs/GCA/009/914/755/GCA_009914755.4/hub.txt;genome=GCA_009914755.4;chrom=CP068254.1;start=5;end=10'
Response:

{ "downloadTime": "2022:12:23T22:41:22Z", "downloadTimeStamp": 1671835282, "hubUrl": "https:\/\/hgdownload.soe.ucsc.edu\/hubs\/GCA\/009\/914\/755\/GCA_009914755.4\/hub.txt", "genome": "GCA_009914755.4", "chrom": "CP068254.1", "start": 5, "end": 10, "dna": "TGTAG"} %
Question: Dear all,

Now we used ANNOVAR to annotate our SNP detected. But the result of ANNOVAR is the single-letter Amino acid abbreviation. It's like this:

The results I want is like this:

So does anyone knows that which annotation software can feed back this 3-letter amino acid abbreviation? Thanks in advance.

Best wishes,

Answer: Hello,

assuming that the annovar output file is a tab separated file, here is a `python` script. Save it as `one2three.py` and run it like this:

$ python one2three.py annovar.csv > annovar_fixed.csv

https://gist.github.com/finswimmer/76bdc20b35ffe5d41bf6a52ff21e1d4b

fin swimmer
Question: Hi all, During Differential gene expression analysis of RNASeq data (DESEq2 or Cufdiff) which is best method to filter differentially expressed genes? Should I go with all the genes having adjusted P value < 0.05 or should I filter them based on a log2 Fold change cut-off?

Thank you
Question: Hello

I have run snp2hla in my dataset and got the imputed results. But I am unable to interpret it. The softwares website also does not have sufficient information. Could anyone help me on this. Thanks in advance.

Answer: The output has P or A against the haplotype which indicates either presence or absence of haplotype
Question: I would like to get sample distances between different samples of an RNA-seq experiment. Read that VST and rlog function of DEseq R package were good to make a correction so that standard deviation of expression of a gene across all samples doesn't change with the mean (of expression of that gene across all samples).
My questions are:

1 - Should these corrections be applied after normalising raw counts for sequencing depth (with the DESeq() function) or directly applied on the raw data?

2 - To do a heatmap with a dendrogram representing the distances between samples, is it better to plot in a heatmap the values corrected with VST/rlog or FPKM values?

3 - 'VST' method seems to be better for big sets (n>30). I have 3 samples, so that means need to choose 'rlog' instead?

4 - In both methods we can set parameter 'blind'. Should I set it to 'TRUE' or 'FALSE' in which situations?

Regards.

Answer: > 1 - Should these corrections be applied after normalising raw counts
> for sequencing depth (with the DESeq() function) or directly applied
> on the raw data?

These should only be applied to the normalised expression levels ('counts'), as per the DESeq2 vignette.

> 2 - To do a heatmap with a dendrogram representing the distances
> between samples, is it better to plot in a heatmap the values
> corrected with VST/rlog or FPKM values?

Don't use FPKM values - the method of normalisation that produces FPKM expression levels should no longer be used for multi-sample studies. Instead, use either the VST- or rlog-transformed counts. Please see the answer that I gave earlier today: https://www.biostars.org/p/328304/#328324

**Edit April 21, 2020:** if you must use FPKM in a heatmap or for any downstream application, I would transform them to z scale via *zFPKM* package

> 3 - 'VST' method seems to be better for big sets (n>30). I have 3
> samples, so that means need to choose 'rlog' instead?

You can justify the use of either. rlog is not recommended for large datasets because it can take a very long time. I tend to check both, where possible, and find that results don't largely change between both of these methods (provided that there are no outliers in your dataset).

> 4 - In both methods we can set parameter 'blind'. Should I set it to
> 'TRUE' or 'FALSE' in which situations?

Please read: http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#blind-dispersion-estimation

Kevin
Question: Hi, in GSE131592, each file corresponds to the counting matrix of one sample. In each file there are three columns. Does anyone know how to get the count of this sample?

The picture shows GSM3790428 as an example. No reply from the dataset author yet.

![example][1]

[1]: /media/images/c34af3bb-79cc-4ea1-94e8-2661c114

Answer: That is htseq count data. The three columns represent the number of reads assigned to that gene assuming the library is 1) unstranded 2) stranded in the forward direction, 3) stranded in the reverse direction.
Question: I am seeing some tools ask for these two types of inversions (head to head or tail to tail), but most callers just produce one breakpoint for each inversion without specifying if it is head-to-head or tail-to-tail. but some callers don't make that distinction and i'm not sure how to extract them. for example, in one delly output i have the following variant:

```
15	48518173	INV00004531	T	.	PASS	PRECISE;SVTYPE=INV;SVMETHOD=EMBL.DELLYv0.8.2;CHR2=15;END=50467460;PE=8;MAPQ=60;CT=3to3;CIPOS=-4,4;CIEND=-4,4;SRMAPQ=60;INSLEN=0;HOMLEN=3;SR=7;SRQ=1;CONSENSUS=TATTGCCACAAATGGATGTGTCAGAGGAGGTAATTAAACTTGTGACCCGCACCATCATTTTTGGGTTTTTGTTCATGAAATCCTTGCATAAGCCAATGTCTA;CE=1.96136;RDRATIO=0.710181;SOMATIC	GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV	0/1:-20.9631,0,-99.9607:10000:PASS:274972:389731:252274:1:47:8:32:10	0/0:0,-12.3096,-141.967:123:PASS:275394:551246:254224:2:61:0:41:0
```
The connection type is 3to3 in this record. From what I assume, since this is an inversion, there should be a second breakpoint that is 5to5. Why isn't there a balancing record for this record or any record that i can find? Perhaps the second breakpoint is directly implied, but in that case, why does delly choose to produce both 5to5 and 3to3 connection types (their totals are in approximately the same proportion), instead of only 5to5 or only 3to3? Maybe I'm not understanding the difference of these types.

Answer: Only simple inversions have both types of paired-end clusters. In cancer you often find complex SVs and in that case these SV types just tell you the type of genomic adjacencies, as described [here][1]

[1]: https://github.com/tobiasrausch/wally#paired-end-view
Question: I am working with fairly low coverage GBS data (average <11 read depth), and as such am wondering if it makes sense for me to remove PCR duplicates from my data, as it seems that these are just adding extra depth. Does anyone know if there is a general rule that should be followed in this situation?

Thanks!

edit: I forgot to mention that I have paired end reads, and I'm not sure if this will change the answer to my question.

Answer: Actually, I think that GBS is one of the very few applications in which you can avoid removing PCR duplicates, because the space that you are sequencing is usually small enough to guarantee that you will find perfect duplicates by chance alone.
However, this depends on several properties. For example, having paired end reads and such a low coverage. you should have few duplicates. If you have a lot, then you have a problem, and your reads represent more the PCR artifacts than the distribution of the sample.
I also suggest that you refer to some of the several software packages develpoed for working on GBS data, such as STACKS: https://github.com/enormandeau/stacks_workflow

EDIT: As Eric correctly pointed out, the correct link to STACKS is this: http://catchenlab.life.illinois.edu/stacks/
I apologize for the mistake.
Question: Title says it all - just looking to set a few of these up for various reasons.

For instance, I may teach a bioinformatics course coming up. I would love to have them all download their data from GEO or somewhere, and there pipeline in a docker container, and then have them independently verify (or disconfirm) a published report.

Thus, if there is a vignette that does the leg-work of identifying the experiment etc. already written that would be best, but honestly just dockerized pipelines (of good quality) would be close - I could make the vignette from there.

Thanks v much!

Answer: You're probably looking for something like https://nf-co.re/ for full, high-quality pipelines.

They are still not trivial to run however.
Question: Hello! I hope somebody can provide some insight into TopHat and Bowtie2 behavior. I have a data set that I aligned to a reference using Bowtie2 and Tophat, both running default settings. When I used Bowtie, the alignment rates were >90% for all samples. With Tophat, however, the rates were all around 45-55%. Since Tophat uses Bowtie, I'm not sure why the resulting alignment rates are so much lower. My best guess is that when Tophat calls Bowtie, it uses different settings than the default settings for a user using Bowtie directly, but I haven't been able to figure this out from the respective manuals. Could anyone explain why this might be? I'm curious because there's apparently something about my reads that is very sensitive to differences in aligners used, and I want to know what it is.

I'm not sure what additional information is needed to answer this question but I'm ready to provide it. Thank you!

ETA: I'm using Tophat2 and Bowtie2, latest versions of both

Answer: Hi everyone, I have an answer to this after consulting with one of the developers. Posting in case anyone here is interested or comes across this when googling.

The difference between how Bowtie2 runs by default and how Tophat2 calls Bowtie2 is in the minimum alignment score required to consider an alignment valid. In Bowtie, this minimum is set to -0.6 - 0.6 * read length. Tophat runs Bowtie with the alignment score minimum set to a constant of -14. So Tophat has a much more stringent criterion for which alignments it considers valid. (When Hisat2 calls Bowtie2, the default minimum threshold is set to 0 - 0.2 * read length. With my samples, this resulted in alignment rates around 75% - right in between the Bowtie and Tophat results.)

The developer I communicated with was surprised that this parameter made such a drastic difference in alignment results, so it's worth paying attention to. This is likely because I'm dealing with diverse wild populations that might have significantly diverged from the reference genome.
Question: Is someone familiar with the internals of CellRanger in terms of how it compares the cellular barcodes obtained from a cDNA library vs the one obtained from a feature barcode library? I am in the situation that I inherited a project where the preprocessing was done with CellRanger, outputting several thousands of "good" cells (so with matched CBs for both cDNA and feature barcode libraries), but I am completely unable to replicate this with indepdendent methods. I did a preprocessing with [Alevin](https://combine-lab.github.io/alevin-tutorial/2020/alevin-features/). This yielded literally no overlap between detected barcodes for cDNA and feature barcode libraries. I then tried aligning the CBs (the first 16bp of R1 of the feature barcode libraries) to the CBs that CellRanger returned resulting in < 1% mapping rate with bowtie2 even in lenient `--very-fast` mode. I verfied with the same strategy that both the cDNA and feature barcode CBs align with > 95% to the 10X-provided 3M-whitelist, so that is not the issue, these are 10X libraries and all other QC is good. Still, no overlap between cDNA and feature barcodes.

Any comments on this? @rob, hope you see this, is there anything fundamentally different towards how Alevin and CR treat/identify CBs?

Answer: Hi @atpoint ,

It's a common problem with feature barcoding pipelines, please checkout the following discussion on the salmon repo.
I'm pretty sure that'll solve the problem, if not please feel free to reach out.

https://github.com/COMBINE-lab/salmon/discussions/576
Question: Hi,

I have 10 fasta files (each file with 20 gene sequences from each of the 10 samples). I would like to create 20 files, specific to each gene from 10 samples.
I proceeded as follows to extract genes with the file_name in header:

pyfasta extract --header --fasta test.fasta gene_name1 | awk '/^>/ {$0=$0 "_sample1"}1' > gene_name1.fasta

Output:

>gene_name1_sample1
ATGC

I am successful in creating multiple gene fasta files for each gene from each sample (a part from loop):

pyfasta extract --header --fasta $sample.fasta gene_name1 >> gene_name1.fasta
pyfasta extract --header --fasta $sample.fasta gene_name2 >> gene_name2.fasta

But, I am unable to add file_name to the header of files in loop (but can do for 1 file as mentioned in the beginning).

Kindly guide.

Thanks.

Answer: The question is pretty generalized with the use of terms. From what I understand, this should work:

If each sample file are not single line fasta, [linearize][1], on commandline:

for file in *.fasta; do awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);} END {printf("\n");}' < $file > "`basename $file .fasta`_single-line.fasta"; done

Then run this python script:

#!/usr/bin/env python

import glob
from collections import defaultdict

genes = defaultdict(list)

for file in glob.glob('*_single-line.fasta'):
with open(file, 'r') as f:
for line in f:
if line.startswith(">"):
# Gene ID : (sample name, sequence)
genes[line.strip().split('>')[1]].append((file.split('.fasta')[0], next(f).strip()))

for g in genes:
# open fasta per gene name
with open(g + '.fasta', 'w') as out:
# for each sample name and sequence
for n in range(len(genes[g])):
# print to file, >gene_sample, sequence
out.write('>' + g + '_' + genes[g][n][0] + '\n' + genes[g][n][1])

[1]: https://www.biostars.org/p/9262/
Question: Hello folks!

I want to mask a genome with a particular repeat library using RepeatMasker.

Then I want to cross the coordinates of the repeats with those of gene annotations to find overlaps between them and study associations and stuff.

I'm only starting to consider feasible ways to do that so any input would be great.

Thanks!

Answer: The answer by Jon is a good one to get overlaps, and I can tell you how to make the GFF. In the 'util' directory of the RepeatMasker distribution there is a script called 'rmOutToGFF3.pl' that will do the conversion. The usage is pretty simple since it writes to stdout:

perl rmOutToGFF3.pl my_repeatmasker.out > my_repeatmasker.gff
Question: Hi,
I am analyzing genomic (RNA-seq) data from Patient Derived Xenograft tumor samples where cancer patient tumors are transplanted and grown in a mouse, harvested, and then extracted DNA and RNA is sequenced. I have never done this before and wondering if human or mouse reference genome should be used. I would guess the human reference genome should be used for alignment but is it possible that there would be some mixed in mouse cells?
Thanks,
- Pankaj

Answer: Hi,
you would use human reference for alignment of course.

There are some methods to remove possible contaminated reads originating from mouse cells. You can use uniquely mapped reads with high mapping quality. There is also a [tool][1] available to separate reads especially from xenograft samples.

[1]: http://bioinformatics.oxfordjournals.org/content/28/12/i172.full
Question: I have chromosome start and end coordinates for a given chromosome. I need to translate these to locus - something like this

"11q1.4-q2.1", meaning it is on the long arm of chromosome 11, somewhere in the range from sub-band 4 of band 1, and sub-band 1 of band 2.
from wikipedia http://en.wikipedia.org/wiki/Locus_%28genetics%29

Answer: You can download the coordinates of the cytobands from the UCSC [here](http://hgdownload.cse.ucsc.edu/goldenPath/hg18/database/cytoBand.txt.gz).

```
curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg18/database/cytoBand.txt.gz" | gunzip -c
chr1 0 2300000 p36.33 gneg
chr1 2300000 5300000 p36.32 gpos25
chr1 5300000 7100000 p36.31 gneg
chr1 7100000 9200000 p36.23 gpos25
chr1 9200000 12600000 p36.22 gneg
chr1 12600000 16100000 p36.21 gpos50
chr1 16100000 20300000 p36.13 gneg
chr1 20300000 23800000 p36.12 gpos25
chr1 23800000 27800000 p36.11 gneg
chr1 27800000 30000000 p35.3 gpos25
(...)
```
Question: Hello everyone,

I would like to have per-base coverage of a bam file of an RNAseq experiment, but only on exons.

I know how to do that on all genome with bedtools:

*bedtools genomecov -d -ibam aligned_reads.bam -g genome.fa > cov.txt*

And I have a gff file giving me the position of exons on my species.

Does anyone know a way to have the same result, but only on given exon positions?

Thanks,

Guillaume

Answer: From the gff extract those lines of exons only. something like below

awkexons.gff

Now use bedtools *coverage* with the above gff and your bam file with -d and other options that suits your purpose
Question: Hello Biostars,

I am very new to genomics and have been given scripts from other bioinformaticians to learn from. Within these scripts they used specific .bed files to analyse a panel of genes and perform annotation.

I understand there are methods of turning a .fasta file or a .bam file into a .bed file. However, I do not understand how to extract information from specific genes e.g. how to create a .bed file which maps all the collagen genes or how to create a .bed file with only exons.

Does anyone know the process of performing such analysis, or know of any databases where .bed files may be stored?

Many thanks,
Krutik

Answer: After [installing BEDOPS][1], here's a way to get a BED file of genes from a central reference like Gencode:

```
wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
| gunzip --stdout - \
| awk '$3 == "gene"' - \
| convert2bed -i gff --attribute-key="gene_name" - \
> genes.bed
```

Likewise, for exons:

```
wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
| gunzip --stdout - \
| awk '$3 == "exon"' - \
| convert2bed -i gff --attribute-key="gene_name" - \
> exons.bed
```

Via: https://bioinformatics.stackexchange.com/questions/895/how-to-obtain-bed-file-with-coordinates-of-all-genes

(Using the `--attribute-key="gene_name"` option with `convert2bed` will bring in HGNC symbol names, which in some contexts can be more commonly used (and useful) for gene names than Ensembl IDs.)

If you want specific genes (say you have a list of gene names or symbols in a file called `genes_of_interest.txt`, you can use `grep`:

```
grep -wfF genes_of_interest.txt genes.bed > genes_of_interest.bed
```

To do mapping of genes to a set of reads of interest, you can use BEDOPS `bam2bed` to convert BAM and `bedmap` to map:

```
bam2bed < reads.bam > reads.bed
bedmap --echo --echo-map-id-uniq reads.bed genes_of_interest.bed > answer.bed
```

[1]: https://bedops.readthedocs.io/en/latest/content/installation.html