cannin/biostars_qa huggingface dataset:
Question: We'd like to perform trio analysis of mice to detect de novo mutations (via WGS and variant calling). The complication: our mouse facility performs breeding in bulk (~15 each males and females). Mom-pup pairs are obviously straightforward. To identify dads, our strategy is to sequence all of the males and infer paternal parentage by comparing dad-pup SNPs. One method is to manually identify SNPs at heterozygous loci for discrimination, but I assume there are tools that can perform this analysis from bulk VCFs. Any suggestions?

Answer: There are several options: I list you a couple that just came to my mind.

The most rigorous one would be to use [Mendel][1]. Method 9 is pedigree selection. Basically you specify all the possible pedigrees and it tells you which one is the most likely. Cons of this approach: 1) I am not sure it is easy to specify you data structure (one trio plus 14 unrelated males) to the software. You might have to do several comparisons between trios where mom-pup are the same and dad is changing. 2) I do not think mendel takes vcf as input.

Another approach would be to use the number of Mendelian inconsistencies to find the most probable father, using for example the [mendelian plugin][2] (I never tried it. For sure GATK has something similar). Basically, you run it on all the possible trios, i.e. for each mum-pup couple you rotate all the possible dads. If the dads are not related to the mom, then one of them (the real one) should show a sensibly lower number of mendelian inconsistencies.

Be advised that either analysis you perform, you should perform it on a very reliable set of SNPs.

[1]: https://www.genetics.ucla.edu/software/Mendel_current_doc.pdf
[2]: https://samtools.github.io/bcftools/howtos/plugin.mendelian.html
Question: Hello everyone,

I am fairly new to bioinformatics. I am currently learning the command line and was writing a very
simple script to align multiple files (Hepatitis C viral DNA sequencing) to reference genome using BWA, and then using samtools to align to the reference.

I seem to be getting a problem in the index step.

I keep getting the same error:

[E::hts_open_format] Failed to open file HCV242.sorted.bam
samtools index: failed to open "HCV242.sorted.bam": No such file or directory

Below is the whole script:

SAMPLE_LIST="HCV176 HCV242"
for SAMPLE_ID in $SAMPLE_LIST

do
#Index the aligner with the reference genome
bwa index -a is HCV1.fasta
# Align reads 1 of samples to ref genome
bwa aln HCV1.fasta $SAMPLE_ID.1.fastq >$SAMPLE_ID.1.sai
# Align reads 2 of samples to ref genome
bwa aln HCV1.fasta $SAMPLE_ID.2.fastq >$SAMPLE_ID.2.sai
# Merge and change into SAM format for both samples
bwa sampe HCV1.fasta $SAMPLE_ID.1.sai $SAMPLE_ID.2.sai $SAMPLE_ID.1.fastq $SAMPLE_ID.2.fastq > $SAMPLE_ID.sam
bwa sampe HCV1.fasta $SAMPLE_ID.1.sai $SAMPLE_ID.2.sai $SAMPLE_ID.1.fastq $SAMPLE_ID.2.fastq > $SAMPLE_ID.sam
# Change the SAM to BAM format for both samples
samtools view -bt HCV1.fasta $SAMPLE_ID.sam>$SAMPLE_ID.bam
# Sort the BAM files
samtools sort $SAMPLE_ID.bam -o $SAMPLE_ID.sorted
# Create the index file for the BAM files
samtools index -b $SAMPLE_ID.sorted.bam
done

I've tested the code and everything works until the index command. All the necessary files are generated as expected up until that step.
Any help much appreciated.

Answer: Are you getting a truncated file name after the `sort` step, like `HCV242.sorted`? In that case you should change `samtools sort $SAMPLE_ID.bam -o $SAMPLE_ID.sorted` to `samtools sort $SAMPLE_ID.bam -o $SAMPLE_ID.sorted.bam` and see if that fixes the issue.
Question: I wish to get all genes present in a particular Reactome Pathway for H. sapiens. Ideally the input would be the names of the pathways or Reactome IDs (e.g.rRNA modification in the nucleus / R-HSA-6790901) and I would like to get back the genes present in that pathway.

Can any bioconductor package do that?

Thanks

Answer: Can use the Reactome Content Service wrapper in R: https://bioconductor.org/packages/devel/bioc/html/ReactomeContentService4R.html

[Event to non-Reactome ids](https://bioconductor.org/packages/devel/bioc/vignettes/ReactomeContentService4R/inst/doc/ReactomeContentService4R.html#Event_to_non-Reactome_ids):

```
library(ReactomeContentService4R)
event2Ids(event.id = "R-HSA-6790901") # can only input id here
```
```
$geneSymbol
[1] "WBSCR22" "MERM1" "HUSSY-03" "PP3381" "BUD23" "TRMT112" "AD-001" "HSPC152"
[9] "HSPC170" "DIMT1" "DIMT1L" "HUSSY-05" "PNO1" "NOP2" "NOL1" "NSUN1"
[17] "RNA45S5" "28S rRNA" "TSR3" "C16orf42" "UND313L" "FBL" "FIB1" "FLRN"
[25] "NOP58" "NOL5" "NOP5" "HSPC120" "SNU13" "NHP2L1" "NOP56" "NOL5A"
[33] "NAT10" "ALP" "KIAA1709" "THUMPD1" "NOC4L" "NOL11" "L14" "UTP18"
[41] "WDR50" "CDABP0061" "CGI-48" "WDR3" "UTP6" "C17orf40" "HCA66" "MHAT"
[49] "WDR36" "TBL3" "SAZD" "PWP2" "PWP2H" "FCF1" "C14orf111" "CGI-35"
[57] "DCAF13" "WDSOF1" "HSPC064" "KRR1" "HRB2" "WDR46" "BING4" "C6orf11"
[65] "FP221" "DHX37" "DDX37" "KIAA1517" "RPS14" "PRO2640" "DDX47" "RPS6"
[73] "OK/SW-cl.2" "EMG1" "C2F" "DDX49" "UTP14C" "KIAA0266" "UTP14A" "SDCCAG16"
[81] "RPS9" "RPS2" "RPS4" "SNORD3A" "RRP9" "RNU3IP2" "U355K" "RCL1"
[89] "RNAC" "RPC2" "RPCL1" "RTC2" "HSPC338" "BMS1" "BMS1L" "KIAA0187"
[97] "PDCD11" "KIAA0185" "NOP14" "C4orf9" "NOL14" "RES4-25" "RRP7A" "CGI-96"
[105] "NOL6" "RRP36" "C6orf153" "HSPC253" "RPS7" "UTP11" "UTP11L" "CGI-94"
[113] "HDCMB12P" "MPHOSPH10" "MPP10" "IMP4" "BXDC4" "IMP3" "C15orf12" "MRPS4"
[121] "UTP3" "CRLZ1" "SAS10" "DDX52" "ROK1" "HUSSY-19" "UTP4" "CIRH1A"
[129] "cPERP-E" "KIAA1988" "UTP15" "HEATR1" "BAP28" "UTP10" "WDR43" "KIAA0007"
[137] "UTP5" "WDR75" "UTP17" "DIEXF" "C1orf107" "DEF" "UTP20" "DRIM"
[145] "NHP2" "NOLA2" "HSPC286" "GAR1" "NOLA1" "DKC1" "NOLA4" "NOP10"
[153] "NOLA3"

...
```
Question: Hi all,

I want to predict the ethnicity information of individuals by using ML method. HapMap genotype data were treated as training set, and unknown individuals' genotype data as testing set. Because this is my first time to do the ethnicity prediction, anyone can give some detail informations (such as methods, typical workflows) or articles making it easy to follow? Many thanks!!!

Joe
Question: Hi,

I was wondering if there was data available online for the trajectories of the different atoms of a protein during protein folding?

I know that Folding@Home does protein folding simulation work - is their trajectory data available online somewhere?

Thanks

Answer: Have a look at the dynameomics project or at MoDEL.
Question: Do you know any public scientific SQL server ?

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29

ENSEMBL http://uswest.ensembl.org/info/data/mysql.html

GO http://www.geneontology.org/GO.database.shtml#mirrors

(I'll give a +1 to each correct answer)

Answer: Don't forget EnsEMBL's new US East MySQL mirror too:

mysql -h useastdb.ensembl.org -u anonymous -P 5306

See here http://www.ensembl.info/blog/2011/07/06/ensembl-now-has-a-second-public-mysql-server/

You can also run your own mirrors in the cloud http://www.ensembl.info/blog/2011/07/12/run-a-private-ensembl-mysql-in-the-cloud/ ;-)
Question: Dear All,
Am I right to use the following script to transform a fastq file (named test.fastq) to a fasta file? THANKS a lot!

#!/usr/bin/perl
use strict;
use Bio::SeqIO;
my $in=Bio::SeqIO->new(-file=>"test.fastq",-format=>'fastq');
my $out=Bio::SeqIO->new(-file=>">test.fasta",-format=>'fasta');
while(my $seq=$in->next_seq()){
$out->write_seq($seq);
}

Answer: If you want to convert fastq to fasta try:

cat test.fastq | paste - - - - | sed 's/^@/>/g'| cut -f1-2 | tr '\t' '\n'
Question: I am working on human protein protein interaction I need to run SVM I need a negative data set comprising of non interacting protein pairs how to generate these.
I found this database Negatome (http://mips.helmholtz-muenchen.de/proj/ppi/negatome/) but the the problem is it talking about most mammals I cant restrict my search to human.
Is there any other database containing human non interacting protein pairs or is there any other method to do this.
Thanks in advance.

Answer: The Negatome database contain some human proteins which you can extract. In addition, there's this [paper][1] (and data [here][2]). Also if your protein interaction data is from AP-MS, you may want to make use of the [CRAPome database][3].

[1]: http://www.sciencedirect.com/science/article/pii/S1046202312001843?via%3Dihub
[2]: http://www.russelllab.org/negatives/
[3]: http://www.crapome.org
Question: I have a script that returns ensembl API info from the variation module. Currently, using the example code it returns the transcript ID of the variant queried, with this code:

my $gene = $gene_adaptor->fetch_by_transcript_id($tv->transcript->dbID);
return $gene->stable_id if (defined $gene->external_name);

I would like to return the gene name, instead of the ENSG#### ID that this command returns. Is anyone familiar with this tool set, and how I might adjust the code to return the gene name instead of the transcript id?

Answer: Hello,

`$gene->external_name` will usually get you the (official) gene symbol. For human genes this would be the HGNC symbol, e.g. "BRCA2".

`$gene->description` will give you the gene name, e.g. "breast cancer 2, early onset".

Note that Ensembl-specific questions can also send to the Ensembl helpdesk (helpdesk@ensembl.org) or the Ensembl developers mailing list (dev@ensembl.org).To the latter you first have to subscribe, though: http://www.ensembl.org/info/about/contact/mailing.html.

Cheers,
Bert
Bert Overduin, Ph.D.
Vertebrate Genomics Team
EMBL - European Bioinformatics Institute
Wellcome Trust Genome Campus
Hinxton, Cambridge CB10 1SD
United Kingdom
http://www.ebi.ac.uk/~bert
Question: I am tying to retrieve mouse mm10 gene information using biomart library in R, but I don't know how to do that

(The information that I need are `mm10.knownGene.name`, `mm10.knownGene.chrom`, `mm10.knownGene.strand`, `mm10.knownGene.txStart`, `mm10.knownGene.txEnd` and `mm10.kgXref.geneSymbol`)

```r
source("http://bioconductor.org/biocLite.R")
biocLite("biomaRt")

library(biomaRt)

mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
listFilters(mouse)

getBM(attributes=c("ensembl_gene_id", "mgi_symbol"), filters= "mgi_symbol", mart=mouse)
```

Answer: ```r
library(biomaRt)
ensembl <- useMart("ensembl", dataset="mmusculus_gene_ensembl")
annot<-getBM(c("ensembl_gene_id", "mgi_symbol", "chromosome_name", "strand", "start_position", "end_position","gene_biotype"), mart=ensembl)

# For older biomart releases I have used
ensembl <- useMart("ENSEMBL_MART_ENSEMBL", dataset="mmusculus_gene_ensembl", host="jul2012.archive.ensembl.org")

#modify the host depending on when your genomic version was released, or which ensembl archive you want to use
```
Question: Hi there :)

I'm working with D. melanogaster, and trying to align Illumina paired-end reads to reference genome using Hisat2. My ultimate goal is to quantify the detected isoforms. However, I have a problem with the output, as the resulting *.sam file has no alignment inside (has only lines starting with @HD and @SQ with no alignments).

The reference genome I use is: `Drosophila_melanogaster.BDGP6.31.dna.genome.fa` (from Ensembl). The annotation I use is: `Drosophila_melanogaster.BDGP6.84.gtf` (also from Ensembl). For hisat2 manual, I'm using: https://ccb.jhu.edu/software/hisat2/manual.shtml#running-hisat2 .

The sequencing center provided me multiple files (70 x 2) in .fastq format. I've renamed them as: `001_R1.fastq`, `002_R2.fastq`, `002_R1.fastq`, `002_R2.fastq` ... etc. `001_R1.fastq` and `001_R2.fastq` are thus paired. Reads were trimmed with cutadapt.

I first indexed the reference : `hisat2-build Drosophila_melanogaster.BDGP6.31.dna.genome.fa`. This worked fine, I have 8x *.ht2 files in my directory.
Then, I extracted the splice sites from the ref: `extract_splice_sites.py Drosophila_melanogaster.BDGP6.84.gtf`. This also worked fine, I have a *.splices.txt in my directory.

Then, and here comes the tricky part, I'd like to run hisat2 iteratively on my *.fastq , defining them as part of pair 1 (-1 parameter) or pair 2 (-2 parameter). As hisat2 takes input files as comma-delimited (from manual > Command-Line > Usage), I tried to run job like this:

hisat2 -x bt2_index.idx -1 `ls
*_R1* | tr '\n' ','` -2 `ls *_R2* | tr '\n' ','` | samtools view -bS > Dmel_hisat.bam

Anyways, this does not seem to be correclty interpreted by hisat2. I don't have error message, but my *.sam contains no alignment.

So, How do you proceed when having multiple paired *.fastq as input?

Many thanks for your help.

Answer: If you have multiple files from the same sample and plan to process them at once anyway then normally you just concatenate them into a single file. If you have multiple samples then **never align them together**. Almost no aligner supports that.

Whenever you have issues with commands like this, use `echo` and a shell script to print out what the exact final command that would be run is and check and see if that seems reasonable. In this case, each of your lists of files ends with a comma, so perhaps that's causing the problem. You might want to use something like snakemake, where you can more easily create list and merge things together into single strings (at least if you know python).
Question: I am trying to use `samtools view` to filter a bam file. I would ideally like to take a sorted bam with the DT flags set by `picard markduplicates` (e.g. DT:Z:SQ for optical/sequencing and DT:Z:LB for PCR/library) and exclude reads with either DT:SQ only OR DT:SQ & DT:LB. As far as I can tell, `samtools` only supports filtering by ***including*** specific tag:value pairs with the `-d` parameter. Is there any workaround here? Or is the only option to write another intermediate file?

Answer: using samjdk http://lindenb.github.io/jvarkit/SamJdk.html

(not tested)

java -jar dist/samjdk.jar -e 'Object att = record.getAttribute("DT"); return att==null || !(att.equals("SQ") || att.equals("LB"));' input.bam
Question: I am using bedtools subtract with large datasets and it keeps crashing, giving the following error

terminate called after throwing an instance of 'std::bad_alloc'

Is there a way to get over this problem in bedtools?

Alternatively is there any other way to find non-overlapping regions for two bed files?

thanks

Answer: bedtools intersect -v
Question: Dear all,

I installed Bio-Linux 8 (based on Ubuntu 14.04 LTS) under VirtualBox 5.0.16 on a Mac OS X 10.9.5 machine. The default user is called `manager`. The problem is that whenever I need root privileges the system asks for a password for `manager` (which I never set). Using default passwords as `password`, `Password`, `Ubuntu` and similar does not help (not even leaving the field blank, or trying my Mac user password). Also, several procedures I found on the web require root privileges, so it is a vicious circle. Any suggestion?

As Bio-Linux is an exquisitely bioinformatics-devoted Linux distro, I hope this question is pertinent here; however feel free to close if not.

Answer: Just found it by chance! I typed `manager` and it worked!
Question: Hello, I have a file which gives me the position of the altered nucleotide or small insertion or deletion (reference allele and altered allele). Is there any tool which can be used to perform these changes on the reference nucleotide sequence and incorporate the ALT_ALLELE changes ?

CHROM	POS	REF_ALLELE	ALT_ALLELE	TRANSCRIPT_ID	GENE_ID	GENE_NAME	REGION	VARIANT_TYPE	REF_AMINO	ALT_AMINO	AMINO_POS	SIFT_SCORE	SIFT_MEDIAN	NUM_SEQS	dbSNP	SIFT_PREDICTION
SL3.0ch00	723860	A	C	mRNA.Solyc00g005060.1.1	gene.Solyc00g005060.1	Solyc00g005060.1	CDS	NONSYNONYMOUS	W	G	52	0	4.32	1	novel	DELETERIOUS (*WARNING! Low confidence)
SL3.0ch00	723867	A	C	mRNA.Solyc00g005060.1.1	gene.Solyc00g005060.1	Solyc00g005060.1	CDS	SYNONYMOUS	G	G	49	1	4.32	1	novel	TOLERATED
SL3.0ch00	723903	T	C	mRNA.Solyc00g005060.1.1	gene.Solyc00g005060.1	Solyc00g005060.1	CDS	SYNONYMOUS	G	G	37	1	4.32	1	novel	TOLERATED

Answer: Give `bcftools consensus` a try

https://samtools.github.io/bcftools/bcftools.html#consensus

or GATK's `FastaAlternateReferenceMaker` (slower than bcftools consensus)

https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_fasta_FastaAlternateReferenceMaker.php
Question: I accidentally deleted my -contigs file generated by an ABySS assembly, but I still have all of the other intermediate files. How can I remake the -contigs file?

Answer: If you rerun your abyss-pe command with the full name of the deleted contigs file as the last argument, it will do whatever is necessary to rebuild the file, i.e.

$ abyss-pe <... your abyss-pe params for assembly ...> myassembly-contigs.fa

or

$ abyss-pe <... your abyss-pe params for assembly ...> contigs

Replace "myassembly" in the first command with whatever the prefix is for your assembly files.

That's the beauty of 'make'. It will only run the commands that are necessary to build the target you specify. Btw, it is always wise to run `abyss-pe` with the `--dry-run` option first, to see what it *would* do before actually running it, i.e.

$ abyss-pe <... your abyss-pe params for assembly ...> myassembly-contigs.fa --dry-run

Btw, `myassembly-contigs.fa` is just a symlink to `myassembly-6.fa`. So if you still have `myassembly-6.fa`, you can recreate `myassembly-contigs.fa` with:

$ ln -s myassembly-6.fa myassembly-contigs.fa
Question: Dear all, please could you advise :

what is the best strategy to call the somatic variants in the absence of germline sample (when only the tumor sample is available) ? thank you,

-- bogdan

Answer: Are you looking for sensitive identification of all somatic variants, particular actionable variants (eg gene fusion breakpoints), or signs of large-scale events such as chromothripsis?

If you're happy to control your FDR by accepting a low sensitivity and ignoring small to medium sized indelsthen it's viable as I demonstrated [here][1]. If you're interesting in TE reactivation then you're going to need a normal as it's practically impossible to determine if those events are germline or not.

[1]: https://rd.springer.com/chapter/10.1007/978-3-319-42044-8_27
Question: Hi,

I need phyloP scpres for the whole human genome. I was able to download 10M positions from http://genome.ucsc.edu/cgi-bin/hgTables using these instructions:

group: Comparative Genomics

track: Conservation

table: your phyloP table of choice

region: genome

output format: data points

However I need the whole genome information, not just 10M.

The website directs me to the downloads page, but I was lost there unfortunately.

Can anyone please link me to the correct place? (Human Grch37 phyloP values in data point format)

Thanks!!

Answer: A bigwig file for phyloP on hg19 whole genome is here http://hgdownload.cse.ucsc.edu/goldenpath/hg19/phyloP100way/

It's a large 8.9Gb bigwig so I do not necessarily suggest converting this into text format all at once, but you could convert to text using bigWigToBedGraph (or convert small chunks at a time with this tool, which accepts -chrom, -start, and -end arguments)

http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/bigWigToBedGraph
Question: I am working with a poorly defined genus of fungi with very few publicly available genomes and I am wanting to cluster the all available ones plus my novel ones into putative species. I am unsure of a good workflow for this.

Would it be better to current phylogenetic workflows of multiple sequence alignment and tree construction and then based on some phylogenetic distance threshold cluster the genomes into species?

OR

Would some form of multiple sequence analysis (with Busco datasets for example) or core genome analysis (with Roary) where we take SNP differences between the MLSA sequences or core genomes and then based on a SNP count threshold cluster the assemblies into putative species?

I apologies if this is a naive question but any help would be appreciated.

Answer: I don't think there is a "best" approach in your situation, given

1) the taxonomy of the group is probably still flaky.

2) being draft genomes, they probably are spotty and have sizeable gaps, and maybe even contaminants. Repeat regions are often left out of such assemblies.

Actually, a good workflow would perform different types of analyses, in order to get a more complete and solid overall picture about the genus. Some suggestions:

[ReferenceSeeker][1] will give you the closest species, in terms of kmer (min-hash) distances and average nucleotide identity (ANI). You can easily find if there are some really close genomes with it.

[Mashtree][2] uses the same kmer (min-hash) distances to group genomes in a dendrogram. The authors don't consider this dendrogram a phylogeny, but I think this dendrogram probably reflects the phylogenetic history of the genomes, anyway. Thus, it may complement a more traditional core genome phylogeny. Whole genome alignments probably aren't a good option, unless all genomes are really close (and, I would argue, also really high quality). Mashtree works as a substitute for whole genome alignments.

A core genome phylogeny would complement the above analyses (but I don't think [Roary works][3] for fungal genomes). As there aren't many genomes available, you can get a SNP-based phylogeny, or even a maximum-likelihood or Bayesian phylogeny with all genes concatenated. Running the same version of BUSCO on all genomes would also give a good dataset for a phylogeny, in addition to informing on the overall quality of the genomes.

[1]: https://github.com/oschwengers/referenceseeker
[2]: https://github.com/lskatz/mashtree
[3]: https://github.com/sanger-pathogens/Roary/issues/425
Question: Dear Biostars,

How can I append sequences of two fastq files ?

Suppose we have two fastq files:

**file1.fastq**

@HEADER
CTCAGNTTGG
+
AAAAA#EEEE
@HEADER
GTGAGTTTAG
+
AA

Answer: I got it, simply:

paste -d '\n' file2.fastq file1.fastq | sed -n 'p;n;n;N;s/\n//p' > result.fastq

Thank you h.mon for you help