cannin/biostars_qa huggingface dataset:
Question: If I have a set of nucleotide sequences of varying lengths and I want to find the number of times a certain subsequences occurs for each sequence in this set, what is the best way to normalize so I can compare the number of hits each sequence has? I feel that longer sequences will (maybe) be more likely to contain more hits so I need to take gene length into account. My first instinct is to just divide by gene length (or by kilobases) to get a frequency per x comparison but I don't know if this is meaningful. How is this kind of comparison normally done?

A related question I have is how I can normalize by position. Say I want to see if in general more hits tend to occur at the beginning of the coding sequence or end of the 5'UTR. How can I normalize so each gene is say, of length 1-100 (including UTRs) that way I can plot a distribution of sequence occurrence by (scaled) position? I know how this could be done but I just want to make sure I do it in a way that is generally considered acceptable.

Answer: Rho index is a generic function calculation in bioinformatics. kralin S. report that the abundance shame of a subsequence can be calculated using the Probability of independent evens:

P(ab)=P(a)*P(b) ===> P(ab….n)=P(a)*P(b)*…P(n)

Rho is set as the ratio between join evens `"ab..n"` and individual evens `a,b,…n`
Then `Rho=P(ab..n)/P(a)*P(b)*P(n)`
`P(a)*P(b)*…P(n)` denoted the expected value of the event ab..n

Example : we have the target sequence of size `k=11"atwnnnvcggg"
`P(a)=`the frequency of 'a', the same manner for t,c and g

P(w)=P(a)+P(t),P(v)=1-P(t) and P(n)=1

Expected (P(atwnnnvcggg))=P(a)*P(t)* P(w)*P(n)* P(n)* P(n)*P(v)*P(c)*P(g) *P(g) *P(g)

frequency(atwnnnvcggg)=number of occurrence of atwnnnvcggg /N-k-1

where N denote the whole sequence size,then;

Rho(atwnnnvcggg)=frequency(atwnnnvcggg)/ Expected (P(atwnnnvcggg))

In seqinr package, rho function doesn't take care the IUPAC coding then we have to edit the function to do the work. This action is not recommended because the number of all possible kmers combinations is too big. To this work, we have to calculate separately:
- The number of occurrence of the target sequence: `gregexpr2` biostrongs package
- Their expected value as denoted previously
[See also Computational Genome Analysis, Richard C. Deonier, Simon Tavaré, Michael S. Waterman][1]

[1]: http://
Question: I want to know the difference between zero and '.' for an Exac value/GnomAD MAF value in an annotated VCF ?
What is the difference between MAF is 0 as not found in the population and MAF is '.' as doesn't exist?
Like for example, when the variant ExAC All MAF value is '.' what does this mean? I have lots of variants with this MAF in ExAC and GnomAD value both in VEP and Annovar annotated VCFs.

Thanks

Answer: `0` is zero.

`.` is for null: data is not available/ unknown /missing
Question: I have a FASTQ file and I'm able to run the FASTQC program to analyse the file. but when I use trim_galore, FASTQC (or the FASTQC option in trim_galore) is not working anymore.

$ fastqc ./sub1_val_1.fq.gz

This is the output:

Started analysis of sub1_val_1.fq.gz
Analysis complete for sub1_val_1.fq.gz
Failed to process file sub1_val_1.fq.gz
java.lang.ArrayIndexOutOfBoundsException: -1
at uk.ac.babraham.FastQC.Modules.SequenceLengthDistribution.calculateDistribution(SequenceLengthDistribution.java:100)
at uk.ac.babraham.FastQC.Modules.SequenceLengthDistribution.raisesError(SequenceLengthDistribution.java:184)
at uk.ac.babraham.FastQC.Report.HTMLReportArchive.startDocument(HTMLReportArchive.java:336)
at uk.ac.babraham.FastQC.Report.HTMLReportArchive.(HTMLReportArchive.java:84)
at uk.ac.babraham.FastQC.Analysis.OfflineRunner.analysisComplete(OfflineRunner.java:155)
at uk.ac.babraham.FastQC.Analysis.AnalysisRunner.run(AnalysisRunner.java:110)
at java.lang.Thread.run(Thread.java:695)

Is the Failed to process file an error because the version is not correct between trim_galore and FastQC?

I found [this][1], [but that wasn't that helpful][2].

I'm using FastQC v0.11.5 and trim_galore v0.4.1.

I subsetted a library (reads in paired-end) using this:

seqtk sample -s100 ./SRR2937435_1.fastq.gz 10000 | gzip > sub1.fastq.gz
seqtk sample -s100 ./SRR2937435_2.fastq.gz 10000 | gzip > sub2.fastq.gz

The sub1_val_1.fq.gz file was after passing sub1.fastq.gz into trim_galore. FastQC with sub1.fastq.gz is working.

[1]: http://seqanswers.com/forums/archive/index.php/t-4846.html
[2]: https://bugs.launchpad.net/ubuntu/+source/fastqc/+bug/1443275

Answer: Update you TrimGalore to at least 0.4.2:

> ### 07-09-16: Version 0.4.2 released
> * Replaced `zcat` with `gunzip -c` so that older versions of Mac OSX do not append a .Z to the end of the file and subsequently fail
> because the file is not present. Dah...
Question: Hi all,

I am working with RNA-seq data of different drug treated human cell lines sequenced on Illumina (2X150 bp chemistry). Data generated is 11-14 Gb.

I am using `HISAT2` for alignment on human genome hg38 build downloaded from [Ensembl database][1] using default parameters. To my surprise, the alignment percentage is very low (`~2 %`) for all the samples.

I have the following observations regarding the data quality :

1. Data quality (phred score) is good i.e. in the range of `32-40`!
2. Illumina adapters were already removed using trimmomatic. No other contaminants.
3. The fastqc metric '[Per base sequence content][2]' fails for all samples. More specifically, A,T,G,C % is deviating considerably at the first 10 base positions.

**Questions**

1. What could possibly be wrong? Any other QC stuff I need to re-consider?
2. More importantly, this is cell-line data, so shall I consider mapping on specific cell line reference OR is it okay if I have just the hg38 standard human reference?

[1]: ftp://ftp.ensembl.org/pub/release-90/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz
[2]: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/4%20Per%20Base%20Sequence%20Content.html

Answer: Hi Vijay,

Do you know how the library was prepared? Strong deviation for the first 10 bases in the nucleotide composition sounds very much like an adapter sequence. Some library preparations (eg. the [Clontech SMARTer Pico](http://www.clontech.com/CN/Products/cDNA_Synthesis_and_Library_Construction/Next_Gen_Sequencing_Kits/Total_RNA-Seq/Pico_Input_Total_RNA_Illumina_Sequencing) kit, amongst others) require trimming from the 5' of sequences before alignment. These adapters are different from the usual illumina sequenceing adapters and will be removed by default by most trimming tools.

If you know the library preparation kit, you can look up the kit's documentation and you'll usually find instructions. For example, from the [Pico docs](www.clontech.com/xxclt_ibcGetAttachment.jsp?cItemId=111761):

> IMPORTANT: The first three nucleotides of the first sequencing read (Read 1) are derived from the template-switching
oligo. These three nucleotides must be trimmed prior to mapping.

But, as you already have the nucleotide content plot you can do without this and just remove the first 10 bases from the 5' of each read. Then your alignment scores should hopefully improve..

I hope this helps!

Phil
Question: Dear Biostars,

I have 80,000 people in a multi-sample VCF which is chunked into roughly 1500 files that span the human genome annotated with VEP including the LOFTEE tool.

I would like to count per chunk the number of High Confidence (annotated as "HC") loss of function variants (LoF) there are per person (either her or hom) per chunk and then add these up so I can see how many HC LoF variants there are per person.

I was initially just using grep to pull out "HC", running bcftools view with the "alt" command and then converting these chunked files into plink files to use it's count function but I was wondering if there was a clean way of doing this directly from a VCF without intermediate steps?

Many thanks for your help

(Similar question to: https://www.biostars.org/p/50923/)

Answer: how about :

bcftools concat -a --file-list your.vcf.txt -O u | bcftools view -O u -i '' | bcftools stats -
Question: Hi everyone,

I am working with mm10 data and using the GRCm38 build 75 GTF from Ensembl. As everyone knows you need the `tss_id` and `p_id` to be present for differential isoform expression (by cuffdiff) when using any GTF other than the cufflinks' merged.gtf. I am using the following command to add the `tss_id` and `p_id` to my ensembl gtf:

```
cuffcompare \
-o cuffcmp \
-C -G \
-r Mus_musculus.GRCm38.75.protein_linc.gtf \
-s mm10.fa \
Mus_musculus.GRCm38.75.protein_linc.gtf
```

To check whether I was doing it correctly, I checked the entries for a particular gene in both the input and output gtfs.

The 'gene' entry for Xkr4 in the original GTF looks like this:

chr1 protein_coding gene 3205901 3671498 . - . gene_id "ENSMUSG00000051951"; gene_name "Xkr4"; gene_source "ensembl_havana"; gene_biotype "protein_coding";

And, these are the entries corresponding to the above coordinates in the output GTF `cuffcmp.combined.gtf`:

```
chr1 processed_transcript exon 3205901 3207317 . - . gene_id "XLOC_000653"; transcript_id "TCONS_00002060"; exon_number "1"; gene_name "Xkr4"; oId "ENSMUST00000162897"; nearest_ref "ENSMUST00000162897"; class_code "="; tss_id "TSS1356";
chr1 processed_transcript exon 3213609 3216344 . - . gene_id "XLOC_000653"; transcript_id "TCONS_00002060"; exon_number "2"; gene_name "Xkr4"; oId "ENSMUST00000162897"; nearest_ref "ENSMUST00000162897"; class_code "="; tss_id "TSS1356";
chr1 processed_transcript exon 3206523 3207317 . - . gene_id "XLOC_000653"; transcript_id "TCONS_00002061"; exon_number "1"; gene_name "Xkr4"; oId "ENSMUST00000159265"; nearest_ref "ENSMUST00000159265"; class_code "="; tss_id "TSS1357";
chr1 processed_transcript exon 3213439 3215632 . - . gene_id "XLOC_000653"; transcript_id "TCONS_00002061"; exon_number "2"; gene_name "Xkr4"; oId "ENSMUST00000159265"; nearest_ref "ENSMUST00000159265"; class_code "="; tss_id "TSS1357";
chr1 protein_coding exon 3214482 3216968 . - . gene_id "XLOC_000653"; transcript_id "TCONS_00002062"; exon_number "1"; gene_name "Xkr4"; oId "ENSMUST00000070533"; nearest_ref "ENSMUST00000070533"; class_code "="; tss_id "TSS1358"; p_id "P1235";
chr1 protein_coding exon 3421702 3421901 . - . gene_id "XLOC_000653"; transcript_id "TCONS_00002062"; exon_number "2"; gene_name "Xkr4"; oId "ENSMUST00000070533"; nearest_ref "ENSMUST00000070533"; class_code "="; tss_id "TSS1358"; p_id "P1235";
chr1 protein_coding exon 3670552 3671498 . - . gene_id "XLOC_000653"; transcript_id "TCONS_00002062"; exon_number "3"; gene_name "Xkr4"; oId "ENSMUST00000070533"; nearest_ref "ENSMUST00000070533"; class_code "="; tss_id "TSS1358"; p_id "P1235";
```

In the output, the gene_id field has XLOC ids instead of Ensembl IDs. Can I fix this to have Ensembl IDs instead? Is there a better way to add `tss_id` and `p_id` to your Ensembl GTF?

Answer: I have developed an Rscript, cuffdiff_gtf_attributes, which can provided the additional attributes `p_id` and `tss_id` as required by cuffdiff to perform all the differential splicing/coding/expression contrasts. I have tested it with Ensembl GTF.
Question: Hi, thanks for taking the time. I am learning Plink these days and find one question. I download the example(called hapmap) in the website and find the number of genotype data in the PED is different from the number of SNP in the MAP. I want to know if anyone knows the reason. Thanks a lot.

Answer: How did you count the genotype data in the PED?

The PED has SNPs for columns and individuals for rows (2 columns per SNP! I.e., `A A` instead of `AA`), while the MAP has SNPs as rows. Furthermore, the first few columns in a PED file are not genotypes - they are

```
Family ID
Individual ID
Paternal ID
Maternal ID
Sex (1=male; 2=female; other=unknown)
Phenotype
```

[Source][1]

Therefore,

wc -l your_map_file.map

should be identical to

awk '{print (NF - 6)/2}' your_ped_file.ped | head -1

The first wc command prints the number of lines in the map file (= number of SNPs in map), the second prints the number of columns (NF) minus the 6 additional IDs, divided by 2 since you always have two columns per SNP.

If your numbers are actually different, try re-downloading the files.

[1]: http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml
Question: I am trying to filter GTEx for liver tissue. I was able to extract the sample names from the "SamplesAtrributesDS.txt" file and find the corresponding column number in the "transcript_expected_count.gct" file for each sample. I then want to use awk to extract the counts from the column number, but I always get this error:

awk: fatal: cannot open file `{print $myvar}' for reading (No such file or directory)

I use the following code:

#!/bin/bash
Genes=$(awk -F "\t" '{if ($6=="Liver") print $1}' GTEx_Analysis_v8_Annotations_SampleAttributesDS.txt)
for gene in $Genes
do
colnr=$(head -n 3 GTEx_Analysis_2017-06-05_v8_RSEMv1.3.0_transcript_expected_count.gct | tail -n 1 | cut -f3- | tr '\t' '\n' | cat -n | grep $gene | cut -f 1)
echo $colnr
if [ -z "$colnr" ]
then
echo "\$colnr is empty"
else
echo $colnr
awk -v myvar=$colnr '{print $myvar}' GTEx_Analysis_2017-06-05_v8_RSEMv1.3.0_transcript_expected_count.gct
fi
done

I would be glad if someone can point out what I am doing wrong.

Answer: I would do something like this, in python, to filter for the columns of interest.

#!/usr/bin/env python3

path_ec = "~/Downloads/GTEx_Analysis_2017-06-05_v8_RSEMv1.3.0_transcript_expected_count.gct.gz"
path_sa = "~/Downloads/GTEx_Analysis_v8_Annotations_SampleAttributesDS.txt"

import pandas as pd
from itertools import count

SKIP_LINES_GTEx = 2
SEP = '\t'

df = pd.read_csv(path_sa, sep=SEP)
ids = df.SAMPID[df.SMTS == "Liver"]

all_cols = pd.read_csv(path_ec, sep=SEP, header=SKIP_LINES_GTEx, nrows=0).columns

# Not all SAMPIDs are in the transcript_expected_count file
ids = ids[ids.isin(all_cols)]

# Columns to keep
cols = ['transcript_id', 'gene_id'] + list(ids)

from tempfile import NamedTemporaryFile as TempFile
with TempFile(mode='w') as out_fd:

# OR:
# with open(output_filename, mode='w') as out_fd:

print(F"Writing to temporary file {out_fd.name}")

# Lines per chunk
CHUNKSIZE = 1000

# Read and write chunkwise
for (df, n) in zip(pd.read_csv(path_ec, sep=SEP, header=SKIP_LINES_GTEx, chunksize=CHUNKSIZE, usecols=cols), count()):
df.to_csv(out_fd, header=(n == 0), sep=SEP, mode='a', index=None)

input("Done. Press ENTER")
Question: I have the nucleotide sequence of 59 genes from strain NC_000913.3 that I want to blast to extract from 20 other strains.

To do this, I am running the below code in a loop in the terminal, where $que is one of the 59 genes and ${i} is one of the 20 strains.

blastn -query $que -db ${i} -outfmt '6 sseq' -max_target_seqs 1

This works fine for the most part, but I am getting issues where the stop codons of NC_000913.3 and the other references are not the same - as they can be TAG/TGA/TAA. I tried removing the stop codons from the query genes, but this caused some of the results to be missing a character.

As an example; I have blasted the sequence for gene [cyaA in NC_000913.3][1] with and without the TGA stop codon against [strain AE014075][2]. The sequence without the stop codon also removes the last nucleotides in the alignment Query: G and Sbjct: T (see below). Is there a way to do this and retain the whole sequence apart from the stop codon (or even better; is there a way for me to do the blast but ignore different stop codons?)

End of alignment with stop codon:

Query 2521 CCGCTATTACAGCAATATTTTTCGTGA 2547
` ||||| || ||||| |||||||| |||
Sbjct 4493873 CCGCTGTTGCAGCAGTATTTTTCTTGA 4493899

End of alignment without stop codon (missing TGA but also G/T):

Query 2521 CCGCTATTACAGCAATATTTTTC 2543
||||| || ||||| ||||||||
Sbjct 4493873 CCGCTGTTGCAGCAGTATTTTTC 4493895

[1]: https://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3?from=3991153&to=3993699&report=fasta
[2]: https://www.ncbi.nlm.nih.gov/nuccore/AE014075

Answer: Blast is a local alignment heuristic. If you need global alignment, that's not what you should be using but something like e.g. needle from the EMBOSS package or exonerate with the affine:global option or ggsearch from the FASTA package.
Question: Hi all,

I'm interested in learning more about best practices for merging RNA-seq technical replicates. I've read many Biostars posts on the matter, but I have a somewhat special case.

**Background:**

I sent RNA samples for sequencing, which were split across 4 lanes each as per normal practice. However, the sequencing depth for this run was much lower than expected and the sequencing core re-sequenced the same samples. Thus, I now essentially have 8x technical replicates per sample, 4x lanes from each of two runs. The repeat runs achieved much better depth.

**Questions:**

1. In this case, should the lower depth data be discarded, or can these data still be used in combination with the updated runs?
2. If combined, is there a need to mitigate batch effect? Aside from read depth, the samples have nearly identical statistics with respect to FastQC analysis, % genome alignment individually with STAR, etc.
3. For combining, what stage is most appropriate? Aside from file sizes, is there any difference between merged .fastq files and merged .bam files? What about at the level of raw counts? In the past, I have merged .bam files from different lanes and found that the effect was summing the raw reads per gene between replicates.

Answer: It's very very unlikely that the library prep was redone. You almost certainly should just put everything together; **running the same library on different lanes or different days does not introduce any technical artifacts.**

The only exception would be if the first run had a serious technical problem with the instrument itself causing the reads to be totally unusable, but I doubt they would have given you the reads at all if that were the case.
Question: I was doing some reading and came across the Macmanes article discussing how minimal trimming is better than stringent trimming for mapability in RNA-Seq reads.

For Trimmomatic I used such settings:

'adapters/TruSeq3-PE-2.fa:2:30:10 HEADCROP:13 LEADING:5 TRAILING:5 SLIDINGWINDOW:4:5 MINLEN:25'

Whilst my sequence per bp quality is high, and my per base sequence content also. My per tile quality was flagged red due to 2 bad tiles (image link below). What would one do, I've seen flow cells with major tile problems across the whole cell. Do 2 cells make a big impact on the data, is there anyway to remove these, or would one bother given the other test parameters passed?

Additionally my kmer content is flagged;however I've read this is typical and not a problem?

Thanks.

https://postimg.org/image/fbtwek7wl/

Answer: Don't worry about a couple bade tiles. The fastQC kmer thing always fails in RNAseq, ignore it.
Question: Hello everyone,
I'm using salmon for the alignment of SE sequencing fastq file.
at some point I have this output:

[warning] Fragment GC bias correction is currently *experimental* in single-end libraries. Please use this option with caution.

It's not clear to me what does it mean to use this parameter ( --gcBias ) with caution, what happens if I use it? can I or not?

Thank you, in particular to who made this amazing tool

Answer: There is a lengthy discussion about it at Salmon issue 83:

["--gcBias" option for single-end reads?][1]

[1]: https://github.com/COMBINE-lab/salmon/issues/83

The crux of the problem is (and I quote **rob-p**, Salmon developer) "Given a single-end read, one does not know the length of the fragment from which it originates". If I am not mistaken, the method implemented just uses the mean fragment length to apply the GC-bias correction, and wasn't subjected to extensive validation, thus it is considered experimental.
Question: How can I extract mutations/variants from Pubmed/PMC articles using terms and without need to download these articles?
Any command line tool I can do that?

Answer: Using [EntrezDirect][1]:

$ esearch -db pubmed -query "osteoporosis" | elink -target snp | esummary | xtract -pattern DocumentSummary -element SNP_ID,NAME,SNP_CLASS,CHRPOS, DOCSUM | head -2
1085307068 GORAB delins 1:170539555 HGVS=NC_000001.11:g.170539556_170539557del,NC_000001.10:g.170508697_170508698del,NG_012237.1:g.12435_12436del,NM_152281.2:c.483_484del,NM_001320252.1:c.-58_-57del,NM_001320252.2:c.-58_-57del,NR_027397.1:n.510_511del,NR_027397.2:n.466_467del,NM_001146039.1:c.483_484del,XM_006711628.4:c.-62_-61del,XM_011510150.3:c.-62_-61del,XM_011510149.2:c.432_433del,XM_017002807.1:c.-62_-61del,XM_024450864.1:c.-58_-57del,XM_011510151.1:c.-62_-61del,NP_689494.2:p.Lys162fs,NP_001139511.1:p.Lys162fs,XP_011508451.1:p.Lys145fs|SEQ=[GA/-]|GENE=GORAB:92344
886042603 COL1A1 del 17:50188957 HGVS=NC_000017.11:g.50188957del,NC_000017.10:g.48266318del,NG_007400.1:g.17683del,NM_000088.3:c.2991del,XM_005257058.4:c.2721del,XM_005257058.1:c.2721del,XM_005257059.4:c.2073del,XM_005257059.1:c.2073del,XM_011524341.1:c.2793del,NP_000079.2:p.Gly998fs,XP_005257115.2:p.Gly908fs,XP_005257116.2:p.Gly692fs,XP_011522643.1:p.Gly932fs|SEQ=[A/-]|GENE=COL1A1:1277

Example based on requested term:

$ esearch -db pubmed -query "HRSV" | elink -target snp | esummary | xtract -pattern DocumentSummary -element SNP_ID,NAME,SNP_CLASS,CHRPOS, DOCSUM | head -3
386834236 GAA snv 17:80104542 HGVS=NC_000017.11:g.80104542T>G,NC_000017.10:g.78078341T>G,NG_009822.1:g.7987T>G|SEQ=[T/G]|GENE=GAA:2548
386834235 GAA del 17:80105111 HGVS=NC_000017.11:g.80105111del,NC_000017.10:g.78078910del,NG_009822.1:g.8556del,NM_000152.5:c.525del,NM_000152.4:c.525del,NM_000152.3:c.525del,NM_001079803.3:c.525del,NM_001079803.2:c.525del,NM_001079803.1:c.525del,NM_001079804.3:c.525del,NM_001079804.2:c.525del,NM_001079804.1:c.525del,XM_005257194.4:c.525del,XM_005257194.1:c.525del,XM_005257193.2:c.525del,XM_005257193.1:c.525del,NP_000143.2:p.Glu176fs,NP_001073271.1:p.Glu176fs,NP_001073272.1:p.Glu176fs,XP_005257251.1:p.Glu176fs,XP_005257250.1:p.Glu176fs|SEQ=[T/-]|GENE=GAA:2548
121907943 GAA snv 17:80118271 HGVS=NC_000017.11:g.80118271C>T,NC_000017.10:g.78092070C>T,NG_009822.1:g.21716C>T,NM_000152.5:c.2560C>T,NM_000152.4:c.2560C>T,NM_000152.3:c.2560C>T,NM_001079803.3:c.2560C>T,NM_001079803.2:c.2560C>T,NM_001079803.1:c.2560C>T,NM_001079804.3:c.2560C>T,NM_001079804.2:c.2560C>T,NM_001079804.1:c.2560C>T,XM_005257194.4:c.2560C>T,XM_005257194.1:c.2560C>T,XM_005257193.2:c.2560C>T,XM_005257193.1:c.2560C>T,NP_000143.2:p.Arg854Ter,NP_001073271.1:p.Arg854Ter,NP_001073272.1:p.Arg854Ter,XP_005257251.1:p.Arg854Ter,XP_005257250.1:p.Arg854Ter|SEQ=[C/T]|GENE=GAA:2548

[1]: http://bit.ly/entrez-direct