cannin/biostars_qa huggingface dataset:
Question: where to get the file of hg19 exon, intron, UTR region ?

I read lots of post in biostar. However it's been a long time post.

Answer: Search for gtf/gff files, these should contain the information you are looking for. Make sure the chromosome notation matches with your reference fasta genome.
Question: How to extract allel, Genotype from vcf file using python or other language for 23GB files?
Well, I am able to right script to get allel but for large VCF file its difficult ? what should other possible way to get allel, Genotype information?

Answer: Extracting genotype information using R.

library(vcfR)
vcf <- read.vcfR(vcf_file, verbose = FALSE )
gt <- extract.gt(vcf, element = c('GT'), as.numeric = TRUE)

For python take a look at the following article.

http://alimanfoo.github.io/2017/06/14/read-vcf.html

Genotypes can also be extracted using `SnpSift.jar` in [snpEff][1] using the following command.

java -jar ../snpEff/SnpSift.jar extractFields annotated.vcf CHROM POS REF ALT "GEN[*].GT" > output.tsv

[1]: http://snpeff.sourceforge.net/
Question: Hello everyone,

I guess it will be a very silly mistake, but I am not able to make this work.

I am using the biomaRt package to obtain the chromosome length of different chromosomes of the human genome. The thing is, when retrieving other information such as the ensembl ID, it works well (using other filters). However, with this code, the programme never stops running. Why? What am I doing wrong?

mart_h <- useMart("ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl")
test <- getBM("chromosome_end", filters="chromosome_name", values=c(1:2), mart_h)
test_2 <- getBM("chromosome_end", filters="chromosome_name", values=1, mart_h)

I know there are more efficient ways of obtaining these lengths, but I would prefer using this package (this code it is part of a pipeline).

Thank you very much, I am pretty sure it will be a very silly thing, but I am not able to solve this.

Have a nice day!

Answer: An alternative and simpler way to obtain chromosome lengths is to use the `GenomeInfoDb` Bioconductor package:

library(GenomeInfoDb)
Seqinfo(genome="hg38")
Seqinfo object with 455 sequences (1 circular) from hg38 genome:
seqnames seqlengths isCircular genome
chr1 248956422 FALSE hg38
chr2 242193529 FALSE hg38
chr3 198295559 FALSE hg38
chr4 190214555 FALSE hg38
chr5 181538259 FALSE hg38
... ... ... ...
chrUn_KI270753v1 62944 FALSE hg38
chrUn_KI270754v1 40191 FALSE hg38
chrUn_KI270755v1 36723 FALSE hg38
chrUn_KI270756v1 79590 FALSE hg38
chrUn_KI270757v1 71251 FALSE hg38

The chromosome lengths are stored in the column *seqlengths*. Take a look at `?Seqinfo` for details.
Question: **I can't which of the following files were used to create the GFF3 files?**

ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/dna/Homo_sapiens.GRCh37.dna.toplevel.fa.gz
ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/dna/Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz

or an entirely different one from the ftp:
ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/dna/

The gff3 file I am using is the following:
ftp://ftp.ensembl.org/pub/grch37/current/gff3/homo_sapiens/Homo_sapiens.GRCh37.87.gff3.gz

Answer: Hi O.rka,

The GFF3 file contains all annotated features on all sequences, including chromsomes, regions not assembled into chromosomes and
haplotype/patch regions.

Therefore, the corresponding FASTA files will be one of the toplevel files here: ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/dna/

Toplevel sequences unmasked:
Homo_sapiens.GRCh37.dna.toplevel.fa.gz

Toplevel soft/hard masked sequences:
Homo_sapiens.GRCh37.dna_sm.toplevel.fa.gz
Homo_sapiens.GRCh37.dna_rm.toplevel.fa.gz

There are further details and descriptions of the file names in the README:
ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/dna/README

Best wishes

Ben
Ensembl Helpdesk
Question: Hi,

I am using FastQC and `EstimateLibraryComplexity` utility from Picard Tools to estimate the complexity of my paired-end RNAseq libraries. I have certain questions about the two tools.

First, according to the documentation of `EstimateLibraryComplexity`,

> Attempts to estimate library complexity from sequence of read pairs alone. Does so by sorting all reads by the first N bases (5 by default) of each read and then comparing reads with the first N bases identical to each other for duplicates. Reads are considered to be duplicates if they match each other with no gaps and an overall mismatch rate less than or equal to MAX_DIFF_RATE (0.03 by default).

So it uses a minimum of 5 bp identity to group reads together and using a max difference rate of 0.03 determines whether they are duplicates or not. Does that mean that if my read length is 100bp and the `MAX_DIFF_RATE` is 0.03, then it will assume reads to be duplicated unless they are different by more than 3 bp?

Second, how many bp does FastQC takes into consideration when determining duplicates? Is there an option in FastQC where you can change the minimum number of bases to be compared between reads in order for them to be considered duplicates.

Thanks!

Answer: > Does that mean that if my read length is 100bp and the MAX_DIFF_RATE is 0.03, then it will assume reads to be duplicated unless they are different by more than 3 bp?

Yes.

> Second, how many bp does FastQC takes into consideration when determining duplicates? Is there an option in FastQC where you can change the minimum number of bases to be compared between reads in order for them to be considered duplicates.

From the manual:

> To cut down on the memory requirements for this module **only sequences which first appear in the first 200,000 sequences in each file are analysed**, but this should be enough to get a good impression for the duplication levels in the whole file. Each sequence is tracked to the end of the file to give a representative count of the overall duplication level. To cut down on the amount of information in the final plot any sequences with more than 10 duplicates are placed into grouped bins to give a clear impression of the overall duplication level without having to show each individual duplication value.

> Because the duplication detection requires an exact **sequence match over the whole length of the sequence**, any **reads over 75bp in length are truncated to 50bp** for the purposes of this analysis. Even so, longer reads are more likely to contain sequencing errors which will artificially increase the observed diversity and will tend to underrepresent highly duplicated sequences.
Question: I have multiple genes collected from diabetic patients, I need to find out the pathways related to these genes using KEGG, Reactome and Wikipathways.
I know how to do that but need to find out if the resulted pathways are related to Diabetes or not.
Please help me to do that, is there a tool or R package help me to find out diabetes pathways ?

Answer: The [Molecular Signatures Database](https://www.gsea-msigdb.org/gsea/msigdb) from Broad/UCSD has diabetes specific terms. Just search for 'diabetes' [here](https://www.gsea-msigdb.org/gsea/msigdb/genesets.jsp) to see what is included. You'll notice in the results that it also aggregates terms from other databases such as KEGG, Reactome, and Wikipathway as well.
Question: I've been recommended to use BUSCO in favour of CEGMA for the assessment of some de novo transcriptome assemblies, but I'm seeing inconsistent results with repeated runs of BUSCO and different numbers of threads (1-6). Despite rerunning with the same parameters I'm getting results which vary by upto 2 BUSCOs. I can't see anywhere suggesting that BUSCO is non-deterministic so am wondering if this is a bug.

I'm using the latest version of BUSCO v1.22, hmmer v3.1b2, python v3.4.3 and blast 2.2.29+ against the eukaryotic dataset in transcriptome mode:
python3 BUSCO_v1.22.py -o test -in transcripts.fasta -l /db/busco/eukaryota -m trans -c 2

Grepping for the complete BUSCOs I see this:

> ./run_test_pe1/short_summary_test_pe1:	302	Complete BUSCOs

> ./run_test_pe2/short_summary_test_pe2:	300	Complete BUSCOs

> ./run_test_pe3/short_summary_test_pe3:	302	Complete BUSCOs

> ./run_test_pe4/short_summary_test_pe4:	301	Complete BUSCOs

> ./run_test_pe5/short_summary_test_pe5:	301	Complete BUSCOs

> ./run_test_pe6/short_summary_test_pe6:	302	Complete BUSCOs

Anyone else seen this? Any comments appreciated.

Answer: Fixed it!

There's a bug in the code where BUSCO reads in the filenames in 'translated_proteins/' with os.listdir(), which returns them in an [arbitrary order][1]. Unfortunately, this order matters for BUSCO to define whether a HMMer hit is complete, duplicated or fragmented - I've not fully pinned down why.

Sorting the list of filenames fixes the inconsistency, and shouldn't have any downstream side-effects. Although, I can't be sure.

It's a trivial change and here's the patch for anyone who is interested. I've also sent it to the author.

--- BUSCO_v1.22.py	2016-07-01 16:10:50.000000000 +0100
+++ BUSCO_v1.22c.py	2016-09-20 15:43:53.850117074 +0100
@@ -833,6 +833,7 @@
if mode == 'trans' or mode == 'transcriptome':
print('*** Running HMMER to confirm transcript orthology ***')
files = os.listdir('%stranslated_proteins/' % mainout)
+ files.sort()
if os.path.exists('%shmmer_output' % mainout) == False:
subprocess.call(['mkdir', '%shmmer_output' % mainout])
group = ''; grouplist = []

[1]: https://docs.python.org/2/library/os.html#os.listdir
Question: I am using EdgeR to find DE genes in my data. However, my data is a bit unusual because, in order to reduce necessary sequencing depth, I have been enriching my sequencing libraries for 300 genes of interest before sequencing and then running EdgeR DE analysis on only those genes. Of the 300 genes, 270 are “test” genes that could be DE between the two experimental conditions and 30 are “control” genes that I am reasonably sure should not change. I am concerned that using EdgeR with the default settings is not appropriate for this analysis because it assumes that most of the genes are not DE. I am also concerned that using only a small number of genes is not compatible with using EdgeR. Are these valid concerns? If so, is there anyway to get around them using EdgeR or another program? Would I have to be worried about the same things if I use DESeq2? Any advice would be appreciated! Thanks!

Answer: Regardless of whether you use edgeR or DESeq2 (or something else entirely), your primary difficulty will be in getting the library size normalization done correctly (I assume you hope that many/most/all of the "test" genes are DE and if these happen to have an asymmetric distribution of fold-changes then you can't use them for this step) . The general strategy is as follows:

1. Create a subset of the data containing only the 30 "control" genes.
2. Determine the normalization factor with this (`calcNormFactors()` in edgeR, `estimateSizeFactors()` in DESeq2).
3. Use the factors determined in #2 with the full dataset.

This ends up being very similar to using ERCC spike-ins. So if it's not clear exactly how to do step 3 then you can search this forum and/or the bioconductor support forum for "edgeR ERCC" or "DESeq2 ERCC" and likely find an example with code (for DESeq2 it's just `sizeFactors(full_dataset) <- sizeFactors(control_genes)`, but I'd have to look up the edgeR syntax).
Question: I have a phased genotype VCF file(i.e haplotype in vcf format). I am interested only in the haplotypes that are there in phased vcf file. Please help or suggest some other way to do it.

Answer: Hello siddazra,

does [bcftools convert][1] with the `--hapsample` option do what you want?

Other solution is mention in this [thread][2].

fin swimmer

[1]: http://www.htslib.org/doc/bcftools.html#convert
[2]: https://www.biostars.org/p/292843/
Question: I have miRNA sequences where the 3' adapter is flanked with 12N(random) UMI sequence+RT primer+Universal adapter.

I have tried using regex function in UMI tools but I am not sure if it will work because of the sequence pattern I have.

Please let me know if I have any chance of processing it with UMI tools.

```
Bold - 3' Adapter || Italics - RT primer || In between both of them -UMI Sequence || After RT primer UNiversal adapter.
NGGGGTGTGACAATGGTGTTT**AACTGTAGGCACCATCAAT**ATACCCTATGCT*AGATCGGAAGAGCACACGTCT*GAA
NTCTGACATGTGTGCG**AACTGTAGGCACCATCAAT**ACGACGACCCTA*AGATCGGAAGAGCACACGTCT*GAACTCCA
```

Answer: You could trim the RT-primer and all sequence to the end of the read, using `bbduk.sh literal=AGATCGGAAGAGCACACGTCT ktrim=r k=15 minlen=0` option. That would ensure that your UMI's will now be at the 3-end most end. Then you should be able to use the `--3prime` option from `umi_tools`.

Final solution based on thread below:

bbduk.sh in=your_file.fq.gz out=stdout.fq literal=AGATCGGAAGAGCACACGTCT ktrim=r k=15 minlen=0 | umi_tools extract --extract-method=string --bc-pattern=NNNNNNNNNNNN --3-prime | bbduk.sh in=stdin.fq out=final.fq literal=AACTGTAGGCACCATCAAT ktrim=r k=15 minlen=0
Question: (cross posted on seqanswers: http://seqanswers.com/forums/showthread.php?t=12542 )

Hi all,

I know, by capillary sequencing, that one of my samples contains a mutation at position:120458492.

Some reads were aligned with bwa and I can clearly see my mutation using samtools tview.

120458491 120458501
CCTGCTCTGGGGAGCTATGCCAGGATGGGTGCC
........R........................
.....C..A..................C.....
........C. ......................
............ ...................
........A..... ...............
,,,,,,,,a,,,,,,,,,,,,,,,, ......
,,,,,,,,a,,,,,,,,,,,,,,,,,,, .
........A........................
........A........................
.............................C...
........A........................
.................................
........A........................
.................................
.................................
........A........................
......................C..........
........A........................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
........A........................
.................................
................................
...............................
.....A........................
,,,,,,,,,gg,,,,,,,,,,,t,,,,
.A........................
.........................

however, we I call the mutations using samtools mpileup

${SAMTOOLS}/samtools mpileup -uf ${HG19} sample.bam |\
${SAMTOOLS}/bcftools/bcftools view -bvcg - > snps.bcf
${SAMTOOLS}/bcftools/bcftools view snps.bcf | gzip --best > snps.vcf.gz

But I can't see the mutation in the VCF.

How can I know why mpileup (or bcftools ?) skipped this mutation ?

Thanks !

Pierre

EDIT:

...and here is the output of pileup (not mpileup):

(...)
chr1 120458492 G 26 C.AaaAA.A.A..A.A,A...A,AA^]. JddfhQacfhhgggahehhhhaB[hg
(...)

Answer: Here is the solution , suggested by swbarnes2 on seqanswers:

Try redoing mpileup with -Buf. That's
worked for me. Like you, I observed
that my sanger-verified mutation
looked fine in pileup, but when I
looked at it in the mpileup run
without -B, the Baq calculations
destroyed the quality scores of that
locus, so it wouldn't call a SNP.
Running it with -B should make your
mpileup look like your pileup at that
locus, and then the SNP will have high
enough quality to be counted.

$ gunzip -c snps.vcf.gz | grep chr1 | grep "120458492"
chr1 120458492 . G A 99 . DP=26;AF1=0.5;CI95=0.5,0.5;DP4=10,2,12,2;MQ=60;PV4=1,0.37,1,1 PL:GT:GQ 255,0,255:0/1:99
Question: Hello

I intend to use Infinium/Illumina **450K human methylation** array data in cancer studies. I need to understand if **data for normal sample** is required? I am planning to download data from tcga (where normal data isn't available).

Also another basic doubt is "if the **normal and tumor sample is mixed** when loading onto the slide"?

Answer: Tumor and normal are not hybridized on the same slide. Whether or not you need normal is not a question of bioinformatics, but of biology. In other words, you need to define what questions you are going to address. If those questions include comparison to normals, you will need normals.
Question: I am trying to estimate memory requirements to run the final steps OMA standalone (after the ALL-vs ALL) on a HPC (slurm), without HOG inference.

From the benchmarks in the manual, the suggested formula is 400MB * pow(nr_genomes, 1.4).

This works for the described metazoan dataset (60 metazoas have been successfully computed using 120GB).

The requirements for bacterial genomes are reported to be lower (50GB for 60 genomes).

I adjusted the formula to match that, so it would be ~166MB * pow(nr_genomes, 1.4).

Does this match other peoples experiences working with bacterial or archaeal genomes?

Would leaving out the HOG inference help me to reduce those requirements significantly?

I was hoping to use a dataset of ~400 genomes but may have to rethink that if it will need around 730Gb of memory.

Answer: Hi Andrew,

I don't think that the scaling behavior of OmaStandalone is the same for bacterial genomes than for eukaryotes. The formula you use gives a rough idea from a few datapoints and are also rather conservative. The memory consumption depends a lot on the size of the genomes and how related they are to each other, i.e. the number of homologs and orthologs.

I'm very confident that you could run 400 bacteria with less than 100GB memory. Would be nice if you can post the amount of memory it required in the end once you're done with the computations.

Deactivating the HOG computation will not significantly reduce the amount of memory (at least not for the bottom-up variant of the algorithm).

Best wishes
Adrian
Question: Hi all,

I'm trying to use Salmon in a differential expression analysis, for this i need a transcriptome fasta.

Given the fact that there is only a genome available for my species, i need to generate the transcriptome (the gene annotation is available for my species, in GFF3 format).

Does anyone know a good way to do that?

Thank you in advance for you advices

Answer: `rsem-prepare-reference` from RSEM ( https://www.biostars.org/p/260395/#260409 ), or `gffread` from StringTie - I prefer the later.
Question: Hi Biostars,

How to set the names of two output files (corresponding to each mate) in `randomreads.sh` script of BBMap package?

Help page specifies the `out=filename` option but when doing like that I get only I file containing both mates. Omitting that option produced two files but with random names.

Here is the command I use:

randomreads.sh ref=C_glabrata_CBS138_current_chromosomes.fasta out=reads.fastq length=75 reads=100 coverage=1 replacenoref=t simplenames=t seed=-1 paired=t metagenome=t

Thanks

P.S. How can I tag [Brian Bushnell][1] ?

[1]: https://www.biostars.org/u/14684/

Answer: You can do following with any additional options you need.

randomreads.sh ref=C_glabrata_CBS138_current_chromosomes.fasta out1=read1.fq out2=read2.fq paired=t addpairnum=t

If you use a single output name then the reads will be interleaved. You can separate them into two files using `reformat.sh`.

reformat.sh in=reads.fq out1=R1.fq out2=R2.fq

You had tagged @Brian the right way.
Question: I am working on the Differential Expression Analysis on RNA-Seq Data between tumor sample and normal sample, while my lecturer suggest me better use the raw expression data instead of using FPKM as expression unit , anyone could tell me why ? or what is the advantage for this ?

Answer: If you want further details, you can visit Joshua Starmer youtube collection of videos

In this one, you can see details about [how DESeq2 normalize RNA-Seq data][1] from raw data

And in this one, [how does it edgeR][2]

[1]: https://www.youtube.com/watch?v=UFB993xufUU
[2]: https://www.youtube.com/watch?v=Wdt6jdi-NQo
Question: Hi all,

I am using the DIAMOND aligner and I have been trying to find what the different fields in the DIAMOND blast tabular file (.m8) are. The fields aren't listed on the DIAMOND github page or paper. Does anyone know what the different tabs are?

Thanks!

Answer: They are supposed to be the same as the default tabular BLAST+, output (`-outfmt 6` or `-m 8`):

'qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore'

...
qseqid means Query Seq-id
sseqid means Subject Seq-id
pident means Percentage of identical matches
length means Alignment length
mismatch means Number of mismatches
gapopen means Number of gap openings
qstart means Start of alignment in query
qend means End of alignment in query
sstart means Start of alignment in subject
send means End of alignment in subject
evalue means Expect value
bitscore means Bit score

**EDIT**: edit down descriptors to the relevant ones, and to fix BLAST+ option.
Question: Trying to run assemblathon_stats.pl to calculate many of the basic contig- and scaffold-level statistics but getting error,

![Screenshot below][1]

[1]: https://s22.postimg.org/cktzl1zr5/Untitled.jpg

Answer: solved by installing perl version 5.10.0 This script does not work under the latest Perl version, need to install an older version (Perl 5.10.0).
Question: Dear biostars community,

I'm having some doubts regarding which measure is ideal for reporting expression values and I thought you could help me with your experience.

I've been dealing with RNA-seq data from two projects now (one single end and other project paired-end reads) and I want to choose an expression measurement suitable to compare samples from different experiments.

From literature (I dig a lot into blogs, papers, etc.. ) and essentially I've summed up the following:

1. Both RPKM and FPKM measures shouldn't be used anymore since they contain an essentially arbitrary scaling factor which is dependent on the average effective length of the transcripts in the underlying sample. Not reproducible, not comparable...
2. TPM measure seems to be more appropriate in dealing with this issue since the sums of normalized reads of each sample are the same across all samples, making it "more suitable" to compare samples. However, its calculation (specifically the denominator term) is also sample dependent and this would be the main reason why I shouldn't use it to directly compare expression values between samples.
3. CPM seems to be a less-normalized measure since it takes into account only library size. On the opposite hand, estimated read count don't normalize samples at all, making it useless to my goal (unless I use some between-sample normalization method).

My point is that TPM seems to be the most reliable expression measurement to compare different samples. Still, TMP performs within-sample normalization (although there's a lot of papers comparing samples based on TPM values).

Do you think TPM is suitable to compare between-samples expression values? If not, which method you would recommend? Should I use any between-sample normalization method?

I'm looking forward to hearing your opinion!
Thanks in advance!

Answer: TPM is not suitable for between-sample normalization because it doesn't account for differences in library composition. It is also very dependent on a few highly expressed genes that may not be the same between your samples.

Instead, you could use the [TMM][1]-normalized counts or the [median of ratio][2] normalization used in DESeq2.

[1]: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-3-r25
[2]: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-10-r106