cannin/biostars_qa huggingface dataset:
Question: We are trying to search for variants in dbSNP but running into difficulties where dbSNP uses a human genome build other than hg19. We have a large number of variants we are examining, and this process needs to be fully automated. We have the chromosome, hg19 base position, and allele. Is there a way to determine the genome build used in dbSNP from results, we are using the entrez eutils docsum.

For example, searching '1[Chromosome] AND 69511[Base Position] AND snp[SNP_CLASS]' gives us this docsum for rs75062661.

http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=snp&id=75062661

However, searching '1[Chromosome] AND 906272[Base Position] AND snp[SNP_CLASS]' doesn't give us rs28507236, but using the GRCH38 position '1[Chromosome] AND 970893[Base Position] AND snp[SNP_CLASS]' does

http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=snp&id=28507236

The CHRPOS field in the first link is related to the GRCh37 genome build, while the CHRPOS field in the second link is related to GRCh38. Is there a way to tell which genome build that field is referring to? We need to search by genomic coordinate, however, right now the only way we can think of to search dbSNP by genomic coordinate it to lift over the coordinates to GRCh38 from GRCH37 and try them both. Even then, we wouldn't know for sure which genome build the resulting SNP is part of. Is there a way to search by build in dbSNP? Or even just to determine which genome build the CHRPOS field is referring to in the dbSNP results?

Answer: dbSNP is using GRCh38 positions in both cases. However, in your first example, the coordinates for GRCh37 and GRCh38 happen to be identical. You can see that by looking on [this page][1]. Note that the HGVS notation for both the NT accessions (which represent chromosome 1 versions) are identical.

[1]: http://www.ncbi.nlm.nih.gov/snp/?term=1%5BChromosome%5D%20AND%2069511%5BBase%20Position%5D%20AND%20snp%5BSNP_CLASS%5D
Question: I am trying to create a taxonomic database using GenomicFeatures package. I downloaded GFF3 file from the NCBI.

Codes :

orf <- GenomicFeatures::makeTxDbFromGFF("orf.gff3",format="auto")

I get following output :

Orf

TxDb object:

Db type: TxDb

Supporting package: GenomicFeatures

Data source: mouse.gff3

Organism: NA

Taxonomy ID: NA

miRBase build ID: NA

Genome: NA

transcript_nrow: 0

exon_nrow: 0

cds_nrow: 0

Db created by: GenomicFeatures package from Bioconductor

Creation time: 2019-05-29 22:32:09 -0500 (Wed, 29 May 2019)

GenomicFeatures version at creation time: 1.32.2

RSQLite version at creation time: 2.1.1

DBSCHEMAVERSION: 1.2

Link to the genome : https://www.ncbi.nlm.nih.gov/nuccore/AY386263.1

As you can see that there are no genes in this database. Can anyone help with this please ?

Answer: Hi [lokraj2003][1],

You can add `gene` features to the gff3 file that you downloaded from [https://www.ncbi.nlm.nih.gov/nuccore/AY386263.1][2],
then re-load it again using the same function.

Like this (it can be any method you like to re-format the original gff3, here for example `awk` with focus on creating `gene` lines, adding `Parent` to each `CDS`, and I leave other detail parsing to you):

```
$ cat orf.gff3 \
| awk 'BEGIN{FS=OFS="\t"} $3!="CDS"{print $0} $3=="CDS"{GENE=$0; gsub("\t0\t", "\t\.\t", GENE); gsub("CDS", "gene", GENE); gsub("cds", "gene", GENE); gsub(";product=.*", "", GENE); print GENE; ID=$9; gsub(".*;protein_id=", "", ID); print $0 ";Parent=gene-" ID}' \
> orf_re.gff3

$ head orf_re.gff3
##sequence-region AY386263.1 1 137241
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=10258
AY386263.1	Genbank	region	1	137241	.	+	.	ID=AY386263.1:1..137241;Dbxref=taxon:10258;country=USA: Iowa;gbkey=Src;genome=genomic;isolate=ORFA;isolation-source=nasal secretions of a lamb at the Iowa Ram Test Station during an outbreak in 1982%2C then passaged in ovine fetal turbinate cells;mol_type=genomic DNA;strain=OV-IA82
AY386263.1	Genbank	gene	2409	2858	.	-	.	ID=gene-AAR98099.1;Dbxref=NCBI_GP:AAR98099.1;Name=AAR98099.1;gbkey=gene
AY386263.1	Genbank	CDS	2409	2858	.	-	0	ID=cds-AAR98099.1;Dbxref=NCBI_GP:AAR98099.1;Name=AAR98099.1;gbkey=CDS;product=ORF001 hypothetical protein;protein_id=AAR98099.1;Parent=gene-AAR98099.1
```

And using that gff3 you'll get:

```
> orf <- GenomicFeatures::makeTxDbFromGFF("orf_re.gff3", format = "auto")
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
> orf
TxDb object:
# Db type: TxDb
# Supporting package: GenomicFeatures
# Data source: orf_re.gff3
# Organism: NA
# Taxonomy ID: NA
# miRBase build ID: NA
# Genome: NA
# transcript_nrow: 130
# exon_nrow: 130
# cds_nrow: 130
# Db created by: GenomicFeatures package from Bioconductor
# Creation time: 2019-05-30 16:57:03 +0200 (Thu, 30 May 2019)
# GenomicFeatures version at creation time: 1.34.8
# RSQLite version at creation time: 2.1.1
# DBSCHEMAVERSION: 1.2
```

Hope it helps. :-)

[1]: https://www.biostars.org/u/46248/
[2]: https://www.ncbi.nlm.nih.gov/nuccore/AY386263.1
Question: Hi all,

I have been trying to use Mutect to compare results from Varscan and other tools. To run MuTect, pre-processing from GATK and Picard tools is necessary.

1\. **Mapped reads using BWA.**

2\. **Convert to sorted BAM using PICARD**

```
java -Xmx4g \
-Djava.io.tmpdir=/tmp \
-jar SortSam.jar \
SO=coordinate \
INPUT=Trimmed_ERR361938_trimmed_bwa.sam \
OUTPUT=Test.bam \
VALIDATION_STRINGENCY=LENIENT \
CREATE_INDEX=true
```

3\. **Mark Duplicates using PICARD**

```
java -Xmx4g \
-Djava.io.tmpdir=/tmp \
-jarpicard-tools-1.119/SortSam.jar \
SO=coordinate \
INPUT=Trimmed_ERR361938_trimmed_bwa.sam \
OUTPUT=Test.bam \
VALIDATION_STRINGENCY=LENIENT \
CREATE_INDEX=true
```

4\. **Realign along INDEL using GATK**

```
java -Xmx4g \
-jar GenomeAnalysisTK.jar \
-T RealignerTargetCreator \
-R /steno-internal/chirag/data/indexGenome/hg19/bwa/hg19.fa \
-o input.bam.list \
-I input.marked.bam
```

**NOW I GET ERROR**

```
##### ERROR
##### ERROR MESSAGE: SAM/BAM file input.marked.bam is malformed: SAM file doesn't have any read groups defined in the header. The GATK no longer supports SAM files without read groups
##### ERROR
```
There is this script which should fix this, but I am not sure of some of the parameter used here,

java -jar ~/unixTools/picard-tools-1.119/AddOrReplaceReadGroups.jar

These parameters need to be used

- RGLB=String
- LB=String Read Group Library Required.
- RGPU=String
- PU=String Read Group platform unit (eg. run barcode) Required.
- RGSM=String
- SM=String Read Group sample name Required.

How do I get information on these parameters, as I am analyzing many published reads.

Are there some other ways to fix this step.

Thanks in advance!

Answer: Hi,

SAM/BAM files need a bit of preprocessing before Picard and GATK can work on them. I'm unable to recall all the steps off the top of my head, but this should help you solve the first problem by adding read groups:

https://www.biostars.org/p/47487/

Remember, dummy read group names will suffice to bypass this error.
Question: Hi

I am trying to convert snp coordinates to rsIDs by using the flag `--update-name`, I am not sure if I am using it correctly, my command is:

plink1.9 --bfile files --update-name snps_coordinates_rs_IDs [1] [2] --out Test_update_name

I am using [1] [2] to dedicate columns.
The snp_coordinate file looks like this. The first column has coordinates the second has rsIDs

1:48899:G:A rs974664996
1:114531:G:C rs918310516
1:133360:T:C rs974682326
1:541281:GGGCCAACCACAGGACA:G rs752188509
1:564782:T:C rs879089824
1:636747:C:T rs974667799

My .bim file has coordinates which I want to convert to rsIDs. I am getting this output

> Error: Invalid --update-name column number.

Am I missing column names
Please assist

Answer: 1. Remove the square brackets from your command line. See the "Interpreting our flag usage summaries" section under https://www.cog-genomics.org/plink/1.9/general_usage .
2. Swap the order of the two numbers: "2 1" instead of "1 2". The [--update-name usage summary][1] is "--update-name [new ID col. number] [old ID col.] [skip]"; the new-ID column number comes BEFORE the old-ID column number. Or you could just remove the column numbers entirely, since "2 1" corresponds to the default behavior (this is mentioned in the first paragraph of the --update-name documentation).

[1]: https://www.cog-genomics.org/plink/1.9/data#update_map
Question: Hello,

When I download release 95 repeats soft masked file "ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.toplevel.fa.gz" it is ~1 Gb, however if I decompress it, it becomes 54 Gb.

This is curious as the same soft repeat masked mouse genome decompressed is 2.7 Gb. Any idea why the human genome is so large and if there is any tool to reformat the fasta file into a smaller one?

Thanks, A

Answer: There are huge numbers of haplotypes in the human GRCh38. In the toplevel DNA sequence files these are represented as the whole chromosome, where most of it is Ns and only the haplotype sequence is actual sequence. This means that the compressed files aren't huge, as they just need to encode how many Ns there are, but decompressed are massive with all the Ns represented.
Question: I have a list of some 200 plant common names that I need to find the respective scientific names for. I could just google every one of them but if possible I want to use a software package that can automate the task (for the sake of reproducibility and my own sanity).

I thought of writing a python script that queries google (for spelling correction and context dependent search) and then grab the supplied Wikipedia URL and then try to find the scientific name with 'beautifulsoup' but turns out that slaps you with captcha pretty fast.

Is there a package/tool/method that comes closest to being an industry standard?

Answer: Using [**EntrezDirect**][1]:

$ esearch -db taxonomy -query "thale cress" | esummary | xtract -pattern DocumentSummary -element CommonName,ScientificName
thale cress Arabidopsis thaliana
$ esearch -db taxonomy -query "zebra fish" | esummary | xtract -pattern DocumentSummary -element CommonName,ScientificName
zebra fish Girella zebra
zebrafish Danio rerio

[1]: http://bit.ly/entrez-direct
Question: Hi all,

I have to run the same program for multiple files.

I would like to know if I can do that using a loop?

Many thanks, friends!

Files: SRR10345445_1.fastq SRR10345445_2.fastq SRR10345446_1.fastq SRR10345446_2.fastq
TrimmomaticPE -threads 30 \
SRR10345445_1.fastq SRR10345445_2.fastq \
SRR10345445_1_PE.fastq SRR10345445_1_SR.fastq SRR10345445_2_PE.fastq SRR10345445_2_SR.fastq \
HEADCROP:12 ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:2:keepBothReads \
SLIDINGWINDOW:4:20 LEADING:5 TRAILING:5 MINLEN:40

TrimmomaticPE -threads 30 \
SRR10345446_1.fastq SRR10345446_2.fastq \
SRR10345446_1_PE.fastq SRR10345446_1_SR.fastq SRR10345446_2_PE.fastq SRR10345446_2_SR.fastq \
HEADCROP:12 ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:2:keepBothReads \
SLIDINGWINDOW:4:20 LEADING:5 TRAILING:5 MINLEN:40

Answer: You just need to loop over one filename, pair the second file to it, and then feed these in to your command. There will be more practical examples elsewhere on the forum, but the general form will be this:

for R1 in /path/to/files/*_1.fastq ; do
R2="${R1%_1.fastq}_2.fastq"
TrimmomaticPE -threads 30 "$R1" "$R2" "${R1%.*}_PE.fastq" "${R1%.*}_SR.fastq" "${R2%.*}_PE.fastq" "${R2%.*}_SR.fastq" ...
done

Haven't tested this, so there might be syntax errors.
Question: Hello

I wish to submit my custom gene sequences of different cultivars of oryza sativa. I was using Bankit to submit my sequence taking the help of this tutorial here at [Bankit][1]

Bankit accepts gene start and end locations and also the corresponding feature names of known features.

I want to submit sequences of genes for which I have found variations in the form of SNPs. How do I submit these entire gene sequences using Bankit?

[1]: https://www.ncbi.nlm.nih.gov/WebSub/html/help/feature-table.html

Answer: See the following links:

Complete Genome Submission Guide

https://www.ncbi.nlm.nih.gov/genbank/genomesubmit/

How to: Submit sequence data to NCBI

https://www.ncbi.nlm.nih.gov/guide/howto/submit-sequence-data/

Submitting Sequences using Specific NCBI Submission Tools

https://www.ncbi.nlm.nih.gov/books/NBK53709/

Formatting your Submission

https://www.ncbi.nlm.nih.gov/books/NBK53702/
Question: Hi,

I am trying to visualise my overlapped chip-seq peak regions which I analysed with Homer mergePeaks function. I have got one venn info file and a "result" file. I would like to use that venn info file then visualise it but when I looked for visualisation libraries or programs, I did not find a method which merits my expectations.

the primary problem is my data is big. ( relatively :) ). I have 19 datasets in one conditions group and 9 datasets in healthy one. I have read making venn diagram for more than 3 datasets would not be smart on biostar tread.

I am trying to find overlapped regions of transcriptional factors that why I wanna know which transcriptional factors sites are most common.

I am a python coding and R mediocre.

Please don't link me to https://www.biostars.org/p/77362/ and https://www.biostars.org/p/66091/ threads I have already read them 0192308 times), also if you think homer is not the best tool for finding the overlaps, please feel free to advice others. (Yes, I do know monkseq)

Thank you very much for your help.

Best regards,
Tunc

Answer: Another alternative is a "binary heatmap". It scales better than a Venn diagram, though combinatorial binding of 19 factors is going to be complex no matter what way you look at it.

Here's a working example: https://gist.github.com/daler/07eb1a95f1e4639f22bd
Question: I have a bam file that I would like sorted karyotypically (not lexicographically) but my contigs are not matching the reference file provided by GATK. Getting the reference file that was originally used for alignment and realigning my sample are unavailable options.

My reference uses "1,2,3,...,X,Y,MT" notation but my bam file uses "chr1,chr2,chr3,...chrX,chrY,chrM" notation. Is there a way to remove the chr prefix and change chrM to MT in my bam file? Can I get by with just revising the header only without messing with the reads in the bam file? Thank you for your help!

Answer: If all your bams are this way, it's probably easier to change your reference to match the bams. (That's just changing a few lines of a fasta and rebuilding an index or two).

If you do end up needing to reformat your bams, there is good advice in previous threads:

- https://www.biostars.org/p/119295/
- https://www.biostars.org/p/13462/
Question: Hi,

I have a table like this;

6 6:29002062:rs7755402 0 29002062 G A
6 6:29004091:rs9468471 0 29004091 A G
6 6:29006250:rs9468473 0 29006250 A G
6 6:29006493:rs9461499 0 29006493 C A
6 6:29006844:rs7743837 0 29006844 G A

I want to remove everything before "rs" in the second column. I know I can use egrep,

egrep -o "(rs\S+)" file | cut -d " " -f 2 > newfile

However, then I`m left with only the rs string;

rs7755402
rs9468471
rs9468473
rs9461499
rs7743837
rs6919044
rs41424052
rs6924824
rs6456886
rs6456887

But I actually want the other columns too.
Any help is greatly appreciated!

Answer: sed 's/\t\b.*:rs/\trs/' file > newfile

Explantion:

s/ = substitution
\t\b = look for tab, and then a word boundary
.*:rs/\trs = ./ matches one or more charaters. you match everything till rs and replace it by tab and rs (ie \trs)
Question: Dear all,

I have pair-end RNA-seq data (Illumina) from parasite and I would like to do De-Novo assembly by TRINITY. I have reference genome of my host organism so I can map my data to host and remove from fastq contaminations.

My plan is:

1. Map with bwa/bowtie/novoaling my pair-end FASTQ files to a host reference genome
2. Remove hits from fastq files (cleaning contaminations)
3. For the rest of FASTQ files use TRINITY for De-Novo transcript assembly

My question is:

May I use aligners (bwa etc.) and align raw fastq files to host DNA and then remove contaminants from fastq files? Question is because my data are from RNA-seq project NOT DNA.

How can I remove the sequences from raw fastq files that align to host DNA (cleaning process)?

Or if you have any other advice how to prepare data to TRINITY pipeline I will appreciate it.

Thank you so much for any comment and sharing your experience.

Answer: I agree with Devon

I would do it in following ways:

1. Map fastq files with tophat2
2. Convert unmapped.bam file to fastq (bamTofatsq) and remap with tophat2, this time provide junctions ( with an option `-j` you got from first run, (if replicates then merge the junctions).
3. The unmapped.bam from this run can be converted to fastq.

I think that this is the fastq from which the reads you are looking for.
Question: Hi everyone
recently I download some SRA files (12 files) using this command from NCBI:

prefetch ERR4172779

all of SRA file were dowloaded properly but when I try to convert them to fastq files using this command :

fastq-dump --split-3 ERR4172779

, it shows this:

2020-08-12T11:42:02 fastq-dump.2.9.4 err: file invalid while opening directory within file system module - failed to resolve column 'READ' idx '1'

2020-08-12T11:42:02 fastq-dump.2.9.4 err: file invalid while opening directory within file system module - failed to resolve column 'READ' idx '1'

2020-08-12T11:42:02 fastq-dump.2.9.4 err: file invalid while opening directory within file system module - failed ERR4172779.sra

==================================================

An error occurred during processing.
A report was generated into the file '/home/mohammad/ncbi_error_report.xml'.
If the problem persists, you may consider sending the file
to 'sra-tools@ncbi.nlm.nih.gov' for assistance.

what's the problem??

Answer: Looking at the ENA archive [here](https://www.ebi.ac.uk/ena/browser/view/ERR4172779), it seems fastq files are not available for this accession. However, they provide BAM files.

(As a side note, in my opinion is better to work with ENA than SRA, see also this answer of mine at https://bioinformatics.stackexchange.com/a/14121/339)
Question: I generated TPM counts from fastq data using salmon. This leaves me with the NM_transcript IDs. I would like to generate the gene symbols from the transcript IDs . Biomart does not recognize transcripts, NCBI Datasets produces an error when I run the entire transcriptome. I have been exploring tximport and tximeta, however, I have run into numerous issues particularly with tximeta not detecting my ref file. Any advice would be greatly appreciated.

Update: I have txiimport and tximeta now running, however they create S4 objects and I am unsure how to make these readable

Answer: tximeta was able to compile all my quant.sf files and summarize to gene level
Question: how many read per sample give me good coverage for soybean rna sequensing with 1.1 GB genome length.
I will use Illumina rRNA depleted plant stranded mRNA library prep for cdna liberary and i like to discover rare gene express
the genome of soybean was sequenced.

Answer: Sequencing depth is usually referenced to be the expected mean coverage at all loci over the target sequences, in the case of RNASeq experiements assuming all transcripts have similar level of expression.

For researchers with a fixed budget, often a critical design question is wether to increase the sequencing depth at the cost of reduced sample numbers or to increase the sample size with limited coverage for each sample.

Necessary coverage is determined by the type of study, gene expression level, size of reference genome, published literature, and best practice defined by the scientific community.

C = LN / G
C : coverage
L : read length
N : number of read
G : haploid genome length

On Hiseq 2500 high output run mode, single flow cell for 2X100bp reads : 4 billion paired end read (claimed by illumina so might expect a bit lower), 0.5 Billion per lane.

C = 2 x100 x 500 000 000 / 1 100 000 000 = 91 X

So you should have about 91X of coverager per lane (~3000$) that you now need to divide into conditions and replicates. IF you haven't heard about randomization, replication, blocking and Sampling, now is a good time to do so.
Question: Is there a better way of downloading the human genome reference sequence in fasta format than downloading it from the UCSC site? BWA protocol asks for an index to be created from the human genome reference multi fasta so I want to get this.
Thanks

[Edited for clarification in response to answers and comments:]

Answer: [The version][1] used by the 1000 genomes project is recommended. The mitochondrial genome in the g1k version is the most widely used [rCRS][2]. The chromosomes and contigs are concatenated, so it is less likely to make mistakes (people frequently concatenate all sequences including different haplotypes from the same region).

We have seen a lot of complications caused by different chromosome names (chr1 vs. 1) or different ordering (chr2 before chr10 or after). It is true that which b37 version to use does not matter too much, but converging to something close to a standard would reduce a lot of unnecessary works for everyone.

[1]: ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/
[2]: http://en.wikipedia.org/wiki/Cambridge_Reference_Sequence
Question: Hi all,

I have a fasta file of an assembled genome of a bacterial strain. I would like to identify the species of this strain (or the closest one possible) using the 16s rRNA sequence.

I tried downloading multiple 16s rRNA sequences from the NCBI database and to perform blastn analysis but I don't think I am using the most efficient method.

Do you have any recommendations on what tool to use to identify the species of a particular strain using 16s rRNA ?

Thank you in advance for your help!

Audrey

Answer: Try this:

Download the [bacterial SSU][1] model from RFAM: https://rfam.org/family/RF00177#tabview=tab9
Then use the .cm file with the tool cmsearch of the [Infernal suite][2] to locate the SSU in the genome, extract the genome sequence at the
best scoring coordinates and submit to NCBI Blast or blast against SILVA. You can repeat that with the LSU as well if you like. If you want to blast the rest of the genome to gain further confidence, then you can use the taxon obtained by SSU search to restrict the blast taxon range and thereby speed up the blast search.

Hope this works out for you.

[1]: https://rfam.org/family/RF00177
[2]: http://eddylab.org/infernal/
Question: how to use birch function in R,birch pakage is removed from CRAN it shows errors while installing from archive and i also want to compare the performance of birch and k means clustering.please anyone help me thanks in advance

Answer: Package 'birch' was removed from the CRAN repository.

Formerly available versions can be obtained from [the archive][1].

Archived on 2014-05-30 as long-standing memory access errors leading to segfaults were ignored.

> Therefore, your performance comparison might be hampered by implementation errors. You will most likely not get any meaningful performance data with respect to the algorithm itself, instead it will be influenced heavily by the implementation problems of birch. If you want a fair comparison you should either fix the package (might be hard) or look for a different implementation.

[1]: http://cran.r-project.org/src/contrib/Archive/birch