cannin/biostars_qa huggingface dataset:
Question: In polyploid genomes with limited variation between subgenomes generally read mapping is a challenge leading to frequent read mismapping and hence calling homeologous variants, eg variations that are actually differences between the subgenomes.

When using discoSnp, you can imagine that these loci could also collapse and variants would appear heterozygous among all samples if the locus is non-branching. One could then use read frequency to decide if the variant is homeologous. One would expect to observe a 50/50 distribution of the two alleles in case of a tetraploid or a 25/75 distribution in case of a true variant as the non variant locus would contribute relatively more to the stack.

However when branching becomes more complex such an approach may become difficult and such variants might still end up in the same locus while actually originating from different loci. Are there any strategies that could be applied to recognize these cases and discern between the different loci?

Related to this, when does discoSnp decide that a graph becomes too complex and decides to split it into separate graphs.

Thanks for the reply.

Answer: Hello Yahan,

Thanks for your question here.

Generally speaking, discoSnp confuses true variants with inter-repetition or inter-genome variations.

There are then several ways of sorting out the true and false positives calls.

- One expects that true variants discriminate individuals and thus discriminate read sets. This is why we proposed the *rank* value, associated with each variant. In our experiments, we have shown that variants with rank < 0.2 are likely false positives, while most (>95%) other ones are true positives (see fig7 of the [paper][1]).
- If variants are mapped on a reference genome (for instance using the VCFcreator tool integrated to disco), uniquely mapped variants are likely to be true positives (marked as *PASS*) while other ones (marked as *MULTIPLE*) are more likely due to inter-repetition or inter-genome variations.

Hoping this helps,

Pierre

[1]: https://www.ncbi.nlm.nih.gov/pubmed/25404127
Question: File1

#CHROM POS ID REF_Zv ALT_lm
chr1A 219620 . T A
chr1A 219648 . A G
chr1A 219867 . A G

file2

#CHROM POS ID REF_Zv ALT_RV
chr1A 219457 . C T
chr1A 219670 . A G
chr1A 219867 . A C

File3

#CHROM POS ID REF_Zv ALT_lm ALT_RV
chr1A 219620 . T A NA
chr1A 219648 . A G NA
chr1A 219867 . A G C
chr1A 219457 . C NA T
chr1A 219670 . A NA C

My command is

awk 'FNR==NR{a[$1,$2];next} {if(a[$1,$2]==""){a[$1,$2]=0};print $1,$2,$3,$4,$5, a[$4,$5]} ' file1 file2 > file3

However,
I can not get the file3 which I want.
Could you help me improve the command?
Thanks,
Fuyou

Answer: Using **R** *merge*:

# example files
file1 <- read.table(text = "#CHROM POS ID REF_Zv ALT_lm
chr1A 219620 . T A
chr1A 219648 . A G
chr1A 219867 . A G", header = TRUE, stringsAsFactors = FALSE,
comment.char = "")
file2 <- read.table(text = "#CHROM POS ID REF_Zv ALT_RV
chr1A 219457 . C T
chr1A 219670 . A G
chr1A 219867 . A C", header = TRUE, stringsAsFactors = FALSE,
comment.char = "")

merge(file1, file2, by.x = c("X.CHROM", "POS", "ID", "REF_Zv"), all = TRUE)
# X.CHROM POS ID REF_Zv ALT_lm ALT_RV
# 1 chr1A 219457 . C T
# 2 chr1A 219620 . T A
# 3 chr1A 219648 . A G
# 4 chr1A 219670 . A G
# 5 chr1A 219867 . A G C
Question: Hi,

From where to download the fully sequenced chloroplast genome of rice (oryza sativa japonica)?

Genome sequence fasta file form ensemble plant and ncbi has only 12 chromosomes with no chloroplast DNA?

Thanks,

Waqas.

Answer: Via http://chloroplast.cbio.psu.edu/ --
http://www.ncbi.nlm.nih.gov/nuccore/NC_001320.1?report=Genbank
Question: ![enter image description here][1]

[1]: https://img.techpowerup.org/190813/screenshot-23.png

Hi,
I am doing reciprocal blast to find ortholog. for that i run i create database and run blastn and retrieve hit sequence fasta through awk 'BEGIN { OFS = "\n" } { print ">"$2, $3 }' blast.out > hit.fa
but when i am trying to align hit.fa to reference genome database, i am getting this issue.

I used this command.--

blastn -query ~/Blast/hit.fa -db ~/Blast/omykiss/genome/omykiss -outfmt 6 -out blast_out -num_alignments 1

please help to solve it.

Thanks in advance

Answer: Hyphens (coding for indels relative to another sequence, from BLAST in this case) are not in {A, C, G, T, N}. You should remove them or replace them with Ns, depending on what you're trying to do.

A [quick Google search][1] reveals this was an issue in some other applications, as well.

[1]: https://github.com/wwood/singlem/issues/17
Question: Suppose I have Protein1, Protein2, Protein3 that are expressed by Gene1, Gene2, Gene3. I want to know Gene1, Gene2, Gene3 in the genome sequence are neighbor or not?

Genome: ========Gene1=====Gene2=======Gene3======

Answer: I suggest you to sort your gff file according to scaffold (chromesome) firstly, and coordinate secondly. they sit at the first column and fourth column, repectively.

REF: http://asia.ensembl.org/info/website/upload/gff.html?redirect=no
Question: I have a set of files with genotypic data, each divided into 3 columns of data, including: MARKER, ID, GENOTYPE.

I would like to use AWK (without changing the order/sorting of the files) to perform a VLOOKUP-like command in order to join the data within the files into a single file as follows:

File1:

BIEC2-99962 HOR_233 G_G
BIEC2-9997 HOR_233 A_G
BIEC2-999748 HOR_233 C_C
BIEC2-999848 HOR_233 G_G
BIEC2-99989 HOR_233 A_A

File2:

BIEC2-9997 HOR_250 A_A
BIEC2-999748 HOR_250 C_C
BIEC2-99989 HOR_250 A_C

File3:

BIEC2-9997 HOR_615 A_G
BIEC2-999748 HOR_615 A_C
BIEC2-999848 HOR_615 A_G
BIEC2-99989 HOR_615 A_C

Expected result:

BIEC2-99962 G_G NA NA
BIEC2-9997 A_G A_A A_G
BIEC2-999748 C_C C_C A_C
BIEC2-999848 G_G NA A_G
BIEC2-99989 A_A A_C A_C

I would appreciate any help on this.

Thanks!

Answer: I strongly recommend you against using awk for complex join operations. I would be able to tell you a million stories, about times when I trusted the output of a awk command, only to discover that it didn't run correctly because of some minor mistake or some exception.

This is how I would do it using R and the dplyr package:

```
> library(dplyr)
> bio1 = read.table('bio1')
> bio2 = read.table('bio2')
> bio3 = read.table('bio3')
> bio1 %>%
left_join(
bio2 %>% # join with the second file (only the first and third column)
select(V1, V3),
by='V1') %>%
left_join(
bio3 %>%
select(V1, V3),
by='V1') %>%
mutate_each(
funs(ifelse(is.na(.), 'NA', as.character(.))), starts_with('V3')) # transform all NAs to the string NA
```
Question: Hi,

I'm a medical student with genetics, and I have come to understand some knowledge in bioinformatics would help me getting a research position in the field. I've been browsing the internet for some course, but all the ones I've found so far did not suit my needs. Does anyone know of an online bioinformatics course that:

1. is affordable
2. is open to an undergraduate medical student
3. provides a certificate
4. teaches you the most important things for genetic, genomic and molecular biology research

Answer: You might want to try a few courses that go through the clinical utility of NGS and bioinformatics (although bioinformatics doesn't stop at NGS it is curently one of the main uses)

Genomic Medicine University of Exeter
[https://www.futurelearn.com/courses/diabetes-genomic-medicine][1]

The Genomics Era: the Future of Genetics in Medicine
[https://www.futurelearn.com/courses/the-genomics-era][2]

Clinical Bioinformatics: Unlocking Genomics in Healthcare
[https://www.futurelearn.com/courses/bioinformatics][3]

Cancer in the 21st Century: the Genomic Revolution
[https://www.futurelearn.com/courses/cancer-and-the-genomic-revolution][4]

They only mention bioinformatics briefly, but they do put Genomic technologies and the analysis in the context of what is useful for the clinician.

Probably the next step is to look at variant interpretation and guidelines - which useful if you're heading down the the route of a Clinical Geneticist specialisation (course is not free - although I think it might be distance learning based).

[Post Graduate Certificate in The Interpretation and Clinical Application of Genomic Data (PGCertICAG)][5]

Also most of the content of the courses seems to be currently aimed at either Rare Disease or Cancer genetics and the other major strand - Microbiology / Virology genomics isn't as well covered.

Bioinformatics courses:

For pure bioinformatics courses not necessarily geared towards clinical areas (mainly Coursera ones as mentioned by previous post):

Bioinformatics Specialization
[https://www.coursera.org/specializations/bioinformatics][6]

Bioinformatic Methods I & II
[https://www.coursera.org/learn/bioinformatics-methods-1][7]
[https://www.coursera.org/learn/bioinformatics-methods-2][8]

Genomic Data Science Specialization
[https://www.coursera.org/learn/python-genomics][9]

Also check out courses on Python (or another easy to learn scripting language such as Ruby), R (or another Stats language), possibly SQL and algorithm development.

[1]: https://www.futurelearn.com/courses/diabetes-genomic-medicine
[2]: https://www.futurelearn.com/courses/the-genomics-era
[3]: https://www.futurelearn.com/courses/bioinformatics
[4]: https://www.futurelearn.com/courses/cancer-and-the-genomic-revolution
[5]: http://www.sgul.ac.uk/postgraduate-certificate-in-the-interpretation-and-clinical-application-of-genomic-data
[6]: https://www.coursera.org/specializations/bioinformatics
[7]: https://www.coursera.org/learn/bioinformatics-methods-1
[8]: https://www.coursera.org/learn/bioinformatics-methods-2
[9]: https://www.coursera.org/learn/python-genomics
Question: Hi,

I am making PCA plot in R using `PCAtools` for different group of individuals. The are 65 groups in the sample metadata. At the moment, for smaller groups I simply specify color for each of them in `colkey` (for instance, `"Group 1" = "#004181"`, `"Group 2" = "#0077BA"`, ................, etc). I am not sure how to add colors for more number of groups and manually specifying them would be tedious. Is there a better way to do this?

Plot_Type = "PCA2D"
pdf(paste0(Plot_Type,"_PC1_vs_PC2.pdf"),height = 7,width = 10)
biplot(p,
x = 'PC1', y = 'PC2',
lab = NULL,
colby = 'Group', colkey = c('Group 1' = '#004181', 'Group 2' = '#0077BA', 'Group 3' = '#4E564D', 'Group 4' = '#6BA3D7', 'Group 5' = '#8B63FF'),
legendPosition = 'right', legendLabSize = 13, legendIconSize = 3.0,
subtitle = 'PC1 vs. PC2')
dev.off()

I will appreciate any help. Thank you!

Best Regards,

Toufiq

Answer: Using *rainbow* to create **n** colours, then *setNames* to create a named vector:

# using 5 as a test, change it to as needed, 65?
n = 5

myCols <- setNames(rainbow(n), paste("Group", 1:n))

myCols
# Group 1 Group 2 Group 3 Group 4 Group 5
# "#FF0000" "#CCFF00" "#00FF66" "#0066FF" "#CC00FF"

# then use within your code:
colby = "Group", colkey = myCols,

**Note:**
- No one can tell a difference between 65 colours.
- There are [other packages/functions](https://stackoverflow.com/q/15282580/680068) to create pallets with **n** number of colours, instead of *rainbow*.
Question: I want to annotate my bacterial genomes and metagenomic samples from gut microflora using the Resfams database available from the Dantas lab website http://www.dantaslab.org/resfams but I am not sure how I can apply it to my samples. I am wondering if anyone knows of a script with commands to use the database. I have fastq and fasta files from raw reads and assembled sequences from both whole bacterial genomes as well as metagenomic samples.

Answer: To annotate your assembled contigs with the Resfams models, you'll need to download [HMMER][1], and produce a 6-frame translation of your contigs, perhaps using a tool like [MetaGeneMark][2]. Use hmmsearch on your fasta of translated sequences and it will output motifs that match to the Resfams models.

Since you already have assembled contigs, your commands might look something like:

gmhmmp -m /path/to/MetaGeneMark_v1.mod -A path/to/output/protein.fasta assembled_contig.fasta
cat path/to/output/protein.fasta | grep -v -e "^$" > path/to/clean_protein.fasta
hmmsearch --tblout path/to/output.tblout.scan /path/to/resfams/models.hmm clean_protein.fasta > /dev/null

Note that the most relevant HMMER output is the tblout file. I usually redirect the full output to /dev/null to reduce clutter. You'll then have to parse the tblout file in whatever way is relevant to your interests. The intermediate step is to remove blank lines from the GeneMark output.

[1]: http://hmmer.org/
[2]: http://exon.gatech.edu/genemark/license_download.cgi
Question: Hi all,

I want to remove all mitochondrial reads from my BAM files, but after doing so the number of mapped reads increased. Here's what I did after indexing:

samtools idxstats .bam | cut -f 1 | grep -v chrM | xargs samtools view -b .bam > clean.bam

However comparing with wc -l says:

samtools view -F 0x904 -c clean.bam
19452291
samtools view -F 0x904 -c .bam
18328474

The second number matches the summary output by tophat2 (fr-frstrand)

any suggestions?

Answer: Your approach looks pretty strange to me, and I would instead do the following:

samtools view -h yourbam.bam | grep -v 'chrM' | samtools view -b > yourbamwithoutchrM.bam
Question: Hello I have a lot of sequences in a FASTA file, and I want to extarct a specific sequence knowing the header ID. for example the header of a sequence is:

NODE_19_length_5758_cluster_19_candidate_1

I know that with `grep` I can extract the header, but i want the below sequences to appear on stdout.

How can I do this on bash?

Answer: The proposed solutions are probably all fine but have the limitation that they first have to iteratively find the correct sequence which can take time if the file is (very) large. The by far simplest solution would be to use `samtools faidx` which in a first step indexes your fasta file and then can use this index to retrieve any sequence in basically no time. The index generation takes like 10 seconds for the entire mouse genome and is therefore not really a limitation, and it only has to be done once per fasta file:

samtools faidx your.fasta NODE_19_length_5758_cluster_19_candidate_1

For many sequences manually:

samtools faidx your.fasta seq1 seq2 seq(...N)

or if you want it automated

(upstream cmd creating "\n"-sep list of seq.names) | xargs samtools faidx your.fasta

Please be sure to browse existing threads, this has been asked many times before.

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

https://www.biostars.org/p/329774/
Question: HI,

I have a network in cytoscape and a list of about 40 genes which i want to extract sub-network of these genes from my network, how I can perform this please?

thank you

Answer: You can check from the wiki [link here][1], that should be handy to perform once you have the network set up. There is Network Analyzer setting which can help you to perform the task. Obviously you have answer here but in case you need a manual to see how to do look at the link.

Also @Jean-Karim Heriche made a similar reply few weeks back about selecting the components from network. Take a look at [this thread.][2]

I would request a bit , please do not think that I am being rude. But just that we already have answer for this question , creating same thread makes it redundant. So before posting if you can try and see threads in Biostars with tags like `cytoscape` and `network` to maximize the usage of the blog to one's benefit. It will be handy for any OP posting here and be more targeted in their query. Just a request.

[1]: http://wiki.cytoscape.org/Cytoscape_3/UserManual#Cytoscape_3.2BAC8-UserManual.2BAC8-Network_Analyzer.Subnetwork_Creation
[2]: https://www.biostars.org/p/186773/
Question: Reference

The chromatin accessibility landscape of primary human cancers, M. R. Corces et al., Science 362, eaav1898 (2018). DOI: 10.1126/science.aav1898

This is a recently published paper in Science. Because I am not studying this sort of things, it's hard for me to understand some of those authors' opinions.

![enter image description here][1]

[1]: https://image.ibb.co/m9vLaA/1.png

The Fig. 4A describes the concept the authors used to identify transcription factor occupancy from ATAC-seq data. In the paper, they say **"TF binding to DNA protects the protein-DNA binding site from transposition while the displacement or depletion of one or more nucleosomes creates high DNA accessibility in the immediate flanking sequence"**. Here I think I can understand it correctly: (1) TF binds to this region, and then (2) it blocks Tn5 binding, so (3) consequently there is a deep notch of ATAC-seq signal, and (4) this notch region tells us where a TF binds to.

Later near the end of this article, the authors said:

**(1) For example, if a noncoding somatic mutation causes the generation of a TF binding site, this mutation could lead to an increase in chromatin accessibility in cis and a concomitant increase in the observed frequency of the mutant allele in ATACseq as compared with that in WGS (Fig. 7A). (2) Similarly, a mutation that inactivates a TF binding site can lead to a decrease in chromatin accessibility and a concomitant decrease in the observed frequency of the mutant allele.**

From Fig. 4A we've known that when TF binds to DNA, the ATAC-seq signal plummets. However, the authors now say *the generation of TF binding site increases chromatin accessibility*. **How could this be possible?** Doesn't the generation of TF binding site make chromatin inaccessible?

Answer: The biological concept is that certain TFs like for example the hematopoietic TF PU.1 binds to the chromatin and opens it up (e.g. by recruiting chromatin remodelers) so that other factors can come in and start their actions like initiating or preparing a region for transcriptional processes. But this then involves the opening of an entire region of like 200-500bp. In contrast, and this is what you are mixing up, is that the TF binding itself protects a tiny part of this open region from Tn cutting events, but this is something like 50bp or so. I attached a plot where you see that in condition blue upon treatment, you see an increased footprint of a certain TF indicating increased binding, but this is a very local event (see the scale on the x-axis). In contrast, the increased binding is also associated with higher overall accessability of the vicinity (see that the flanking regions of this 200bp window show increased readcounts on the y-axis). So do not mix up chromatin accessablity induced by TFs with local binding site protection due to the presence of the TF itself. The latter is simply the consequence of the physical presence of a protein denying access to the ATAC-seq Tn.

![enter image description here][1]

[1]: https://image.ibb.co/etbwFA/Screen-Shot-2018-11-06-at-11-22-04.png
Question: Hello everyone !

I'm trying to polish my genome assembly, and for this I need to use Pilon (there other options from PacBio tools and software being unavailable). My problem is that I don't see anywhere some recommendations to install this tool.

I read the several sections form the wiki page : https://github.com/broadinstitute/pilon/wiki , but wasn't able to find it. I found an installation using conda that didn't worked for me.

The steps I accomplished so far, is cloning the pilon deposit from github in my computer. Then, I tried to launch the build.sh script, which gave me the following error :

Loutre:~/pilon$ ./build.sh
./build.sh: 14: ./build.sh: sbt: not found
ln: impossible de créer le lien symbolique «/home/loutre/lib/pilon/pilon-dev.jar»: Aucun fichier ou dossier de ce type

(It says that the file or directory does not exist, which is true in my case).

I think that I missing something. As it's a java source file, I thought about using jar with some other arguments, but there is no single indications in the README.md, and there is so many .scala files in the directory src. I'm kind of lost here.

I saw that some users recommended Pilon to other people, someone also had the same trouble as me ( smkopac) in this post : https://www.biostars.org/p/169796/ when trying to run :

java -Xmx16G -jar lib/htsjdk-1.130.jar pilon --help

Any help or further informations that I may have missed will be really appreciated. I hope that this is not a stupid question.

Cheers,

Roxane

Answer: I allow myself to add my own answer to this question, in case some other people wonder.

So, clone the deposit is probably not the way we should install Pilon.

I've ran successfully Pilon after directly the .jar file available here : https://github.com/broadinstitute/pilon/releases and then typing on a terminal, where your jar file is located :

java -jar pilon-1.21.jar --help
Question: I would like to annotate records in one VCF file (`input.vcf`) with some of the INFO fields of the corresponding records from the database (`db.vcf`), but only if the recorded mutation matches exactly in input and in the database. E. g. let's say I have three very simple VCF files:

`input.vcf`

```
##fileformat=VCFv4.1
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 878638 . G A 100 PASS A=3.0
```

`db1.vcf`

```
##fileformat=VCFv4.1
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 878638 . G A 100 PASS B=4.0
```

`db2.vcf`

```
##fileformat=VCFv4.1
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 878638 . G T 100 PASS B=4.0
```

Note that db1 and db2 describe different SNPs at the same locus; SNP in `db1.vcf` matches with the one in `input.vcf`, but SNP in `db2.vcf` does not. I need a tool that can discern such cases and annotate the input file record with information from database only if the mutations match. Is there a tool to accomplish what I want?

I tried using GATK's VariantAnnotator and vcflib's vcfaddinfo; they unfortunately both ignore information about the mutation and add B=4.0 annotation in both cases.

Just to clarify, this is what I want in the case described:

```
$ some_tool input.vcf db1.vcf # SNP in input and database matches
##fileformat=VCFv4.1
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 878638 . G A 100 PASS A=3.0;B=4.0
```

```
$ some_tool input.vcf db2.vcf # SNP in input and database do not match
##fileformat=VCFv4.1
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 878638 . G A 100 PASS A=3.0
```

Answer: Try bcftools annotate

bgzip -c input.vcf > input.vcf.gz; tabix input.vcf.gz
bgzip -c db.vcf > db.vcf.gz; tabix db.vcf.gz

bcftools annotate -a db.vcf.gz -c CHROM,POS,REF,ALT,INFO/B input.vcf.gz > output.vcf

this will fill in INFO/B from db.vcf.gz when all of CHROM,POS,REF and ALT match.