cannin/biostars_qa huggingface dataset:
Question: A colleague has just found a Human mutated sequence having a Valine (GTG ) as the initiation codon, instead of the very classical Methionine (ATG).

She wonders if this protein could be normally expressed.

I found that Valine could be an "Alternative Initiation Codons" ( see http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi#SG1 ) but I'd like to know if some other ATG in the sequence could be some good candidates for an alternative initiation of the translation.

What is your favorite web server for doing this kind of prediction ?

Thanks,

Pierre

Answer: I cannot offer a tool to use but a person to contact. You should contact Sarah Calvo (of Vamsi Mootha's group at the Broad Institute). She has looked very deeply into 5'-UTRs and upstream open reading frames. In many genes, there are short ORFs that initiate translation from a sequence distinctly upstream from the true Met1 initiation codon. After seeing Sarah speak on more than one occasion, I feel that she would be able to offer some solid advice on the research question at hand.
Question: Hi all,

I'm dealing with a fasta file with spaces at the end of line, which caused the problem. I didn't find a suitable way to remove them. Please kindly tell me the appropriate command for removing them?

Answer: Is that you have an space or a lack of the end of line code?

If your data are tab separated, and you have an space only at the end of the lane, you can do the following

cat file.fasta | tr -d " " > newfile.fasta

But notice that this will get rid of all spaces, including those at the middle of the lane.
Question: I have a list of SRA accession numbers that I download using `sratoolkit`'s `fasterq-dump`. Since I have a number of samples, instead of downloading them serially, I take advantage of array jobs. The script I use is the following:

#!/bin/bash
#$ -cwd
#$ -V
#$ -t 1-44
read -a samples <<< $(cut -d , -f 1 linker.csv | tail -n +2)
false_index=$SGE_TASK_ID
true_index=$((false_index-1))
sample=${samples[$true_index]}
fasterq-dump $sample -O raw_samples >> downloading.log

Short version is that I read the files to download into an array. I specify the amount of threads with `-t` flag. Every task gets id (`$SGE_TASK_ID` variable) from 1 to 44. Based on that id, sample is retrieved via indexing, and then I run the command `fasterq-dump` to download it. This script works wonders.

But I want to implement it in a `snakemake` pipeline, and in order to make that pipeline work, I need to use it's `-cluster` flag. I've managed to write one `snakemake` rule to download files, but only serially. I have no idea how I'd implement parallelization (technically, an array job). The `snakemake ` script is below:

import csv
def read_dictionary():
with open("linker.csv") as csv_file:
csv_reader=csv.reader(csv_file,delimiter=",")
dic = {row[0]:row[2] for row in csv_reader}
return(dic)

SRA_MAPPING = read_dictionary()
SRAFILES = list(SRA_MAPPING.keys())[1:]

rule download_files:
output:
"raw_samples/download.log"
run:
for file in SRAFILES:
#shell("touch raw_samples/{file} >> raw_samples/download.log")
shell("fasterq-dump {file} -O raw_samples >> {output}")

In a nutshell, it reads the samples into global variable `SRAFILES`, and then runs the python code that selects a file from said variable, then runs `fasterq-dump` . How would I implement "parallelization" of one job/rule?

Answer: So `SRAFILES` is a list of SRR IDs, right? If so, I would make the SRR id a wildcard and let snakemake handle the for-loop. E.g.:

```
SRAFILES = ['SRR123', 'SRR456']

rule all:
input:
expand('raw_samples/{srr}.log', srr= SRAFILES),

rule download_files:
output:
xlog= 'raw_samples/{srr}.log',
shell:
r"""
fasterq-dump {wildcards.srr} -O raw_samples > {output.xlog}
"""
```

Then execute with `-j` option set to a reasonable number of jobs to run in parallel:

```
snakemake -j 10 -p --cluster "qsub ..."
```

Basically, you are asking snakemake to produce one log file per SRR id and this log file is produced by rule `download_files`. As a side product, `download_files` will give you the actual fastq files (things could be done differently in this respect but hopefully this will help...)
Question: Does anyone can help explain what it means by the different colors of highlighted reads? Thanks

![Reads on IGV][1]

[1]: /media/images/53d1c555-9bb1-43dc-9eaf-ef6fc524

Answer: Coloring of Mapped Paired Reads in https://software.broadinstitute.org/software/igv/AlignmentData

https://software.broadinstitute.org/software/igv/interpreting_insert_size

https://software.broadinstitute.org/software/igv/interpreting_pair_orientations

https://software.broadinstitute.org/software/igv/interpreting_bisulfite_mode
Question: I read the Refseq documentation in the NCBI handbook but it is still not clear to me.
I'm aware WGS represents all assembled contigs from a sequencing project, and Refseq supposedly has some curation...

Does that mean WGS is more complete than Refseq (even if it includes a bunch of unnannotated features)?

Answer: I'd not think those things are comparable really, as they mean different things. Annotation is only possible when the sequences are available. RefSeq and others provides the annotation of these sequences (e.g. the [Ensembl gene set][1]), whether they are assembled or not (yet). The genomic sequence comes from Whole Genome Sequencing (WGS) experiments and we carry out the annotation of genes, transcripts, genetic variants, regulatory regions, etc.

[1]: http://www.ensembl.org/info/genome/genebuild/genome_annotation.html
Question: Hello,

I have several .fq files containing 5bp inline barcodes at the beginning of each read such as (barcodes are between *) :

@gi|110640213|ref|NC_008253.1|_418_952_1:0:0_1:0:0_0/1
*CCAGG*CAGTGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCATCTGGTAGCGATGAT
+
2222222222222222222222222222222222222222222222222222222222222222222222
@gi|110640213|ref|NC_008253.1|_31_476_0:0:0_0:0:0_1/1
*CAGAT*GGTTGGTGATTTTGGCGGGGGCAGAGAGGACGGTGGCCACCTGCCCCTGCCTGGCATTGCTTTCC
+
2222222222222222222222222222222222222222222222222222222222222222222222
@gi|110640213|ref|NC_008253.1|_210_743_2:0:0_1:1:0_2/1
*CATTA*CCACCACCATCACCATTACCACAGGAAACGGTGCGGGCTGACGCGTACAGGAAACACCGAAAAAA
+
2222222222222222222222222222222222222222222222222222222222222222222222

I would like to modify these sequences in order to have the same for each read (here it would start by AAAAA):

@gi|110640213|ref|NC_008253.1|_418_952_1:0:0_1:0:0_0/1
*AAAAA*CAGTGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCATCTGGTAGCGATGAT
+
2222222222222222222222222222222222222222222222222222222222222222222222
@gi|110640213|ref|NC_008253.1|_31_476_0:0:0_0:0:0_1/1
*AAAAA*GGTTGGTGATTTTGGCGGGGGCAGAGAGGACGGTGGCCACCTGCCCCTGCCTGGCATTGCTTTCC
+
2222222222222222222222222222222222222222222222222222222222222222222222
@gi|110640213|ref|NC_008253.1|_210_743_2:0:0_1:1:0_2/1
*AAAAA*CCACCACCATCACCATTACCACAGGAAACGGTGCGGGCTGACGCGTACAGGAAACACCGAAAAAA
+
2222222222222222222222222222222222222222222222222222222222222222222222

I want to make sure that only the sequence at the beginning of the reads are modified and not throughout the read itself. The barcode sequence might be present within reads and I don't want to modify it.

Do you know any easy way to do this?
Thanks!

Answer: I assume that the * are not part of the sequence and are just there to highlight them :-) Then use awk:

zcat [in fasta file] |awk '{if(NR%4==2){print "AAAAA"substr($0,5)}else{print $0}}' |gzip > [output fasta].gz
Question: Is there any way to detect if your rna-seq data is unstranded or stranded

Answer: A few RNA-Seq QC tools will detect whether a run is strand-specific. For example, the `infer_experiment.py` script in the following claims to do this (never used this myself, so can't vouch for it):

http://rseqc.sourceforge.net/
Question: Hi,
I have low-depth genome sequencing data from several species in a family and I would like to build gene trees as input for analyses with programs like BUCKy and Astral. (I want to estimate a species tree, and assess evidence for introgression between nonsister lineages). I have mapped the reads to the closest outgroup available and produced VCF files for each scaffold.

I am wondering whether "gene trees" need to be real genes, or whether it is acceptable to simply cut the VCF file into small pieces and build "gene" trees of these pieces? It seems like most papers use real genes however I do not have a transcriptome for my organisms unless I used a fairly distant outgroup as a reference to assemble the gene sequences.

Answer: Gene trees don't need to be generated from 'genes' in the sense that they're absolutely restricted to coding loci, but you should be working at a fine enough scale where you'd expect a similar model of sequence evolution and linkage to apply. I've done this estimating gene trees for loci concatenated across chromosomes and by reducing the scale down to exons within a 100kbp window or so, looking for any discordance in support and results once thrown into species tree analyses.

Be careful, though, if you're using model-based phylogenetics to estimate your gene trees - you'll want to include confidently called reference alleles as well to flesh out your alignments and allow for proper estimation of a proportion of invariable sites and/or parameters in the gamma distribution describing among-site rate heterogeneity. Multisample VCFs, once converted to FASTA/Nexus to estimate your gene trees, will only include variable positions and may not be appropriate. You'll also want to be extremely careful with your filtering criteria, especially if you're working with low-coverage data that could introduce biases.

Another point of advice: it may be easier to manipulate FASTA files programmatially as opposed to manipulating the VCFs early on since you ultimately need to work with multispecies alignments per gene/locus/region/scaffold. Check out SeqIO in Biopython and Biostrings in R (part of Bioconductor).
Question: TCGA has recently migrated to the Genomic Data Commons (GDC). Following this migration, many tools convenient for retrieving TCGA data, such as TCGA-Assembler, no longer apply. So, with the new GDC, I'd like to download RNA-Seq data (in bulk) for tumor samples as well as normal control samples. How might I accomplish this?

I know how to download data from the GDC, but I need to know whether a specific RNA-seq data file is coming from a tumor how to obtain its "matched normal tissue" RNA-seq data. I'd like to do this in bulk, for, say, the HCC project.

I know there's a "legacy portal" for the old TCGA data on the GDC website, but I want to use the newest GDC portal.

Thanks in advance.

Answer: Hi everyone,

After looking around, it seems that a few people are having similar problems. So, I contacted GDC and they responded very quickly and helped me find a solution, which I'm posting here in hopes that it well help someone in the future :)

Here's how you do it:

1) Select all the files you want to download, and get the manifest file

2) Use the GDC api (documented here: https://gdc-docs.nci.nih.gov/) as follows: Use the "filters" parameter to get only files with files.file_id matching the UUIDs of the files you want to download (those UUIDs are the first column of the manifest file). Use the "fields" parameter, and set it to "file_id,file_name,cases.submitter_id,cases.samples.sample_type" -- this will get you the file name, the patient (i.e. cases.submitter_id), and what type of sample is it (tumor, normal, etc.)

Here's some quick code (which the person at the GDC help center was kind enough to provide):

{
"filters":{
"op":"in",
"content":{
"field":"files.file_id",
"value":[
"9c5e4668-5ed8-4d7e-972a-6e985b0030df",
"629fb40f-e97d-4bb8-8f59-0a71256211b1"
]
}
},
"format":"TSV",
"fields":"file_id,file_name,cases.submitter_id,cases.samples.sample_type",
"size":"100000"
}

Replace UUIDs with the UUIDs that you're interested in. With some basic shell and JSON, it should be pretty straightforward :)
Question: Is this behavior a result of strand specificity? For example, on the first line would "FUN_003800-T1" be on the negative strand while "FUN_003804-T1" on the positive strand so the `ssend` would be a smaller value than the `sstart`?

![enter image description here][1]

[1]: https://i.imgur.com/Gn5Aukn.png

Answer: Yes, the match is on the negative strand.
Question: Hi,

Any idea if there's a way to get the Instrument type (i.e. NextSeq, HiSeq, etc), from the instrument name field in the Fastq read header? - I've tried looking for a map of sorts in Illumina's documentation, but to no avail.

Thanks,

Answer: I realize that this is a year and a half late, but [here is some code](https://github.com/10XGenomics/supernova/blob/master/tenkit/lib/python/tenkit/illumina_instrument.py#L12-L45) from 10X that incorporates both the flow cell ID and the machine ID to figure out the run type.
Question: Hi all,

The current version of MAF is 2.4. Does MuSiC support this or has any issues with this version, if we make a mutation file that conforms to its specification?

Are all of the MAF columns required? If not, which ones can be excluded and is the ordering of the columns important?

Many thanks in advance.

Answer: Nope, to both questions.

The changes from MAF v2.3 to v2.4 are [detailed here][1]. Of the columns that MuSiC cares about, only `Variant_Classification` has some changes - the values `De_novo_Start_InFrame` and `De_novo_Start_OutOfFrame` are no longer allowed. But MuSiC still treats them as valid values for backward compatibility (see [this line of code][2]).

The order of the first 34 MAF columns is important. But only the following columns are actually used by the tools in MuSiC, while the others can be left blank, or can contain gibberish. `Hugo_Symbol`, `Chromosome`, `Start_position`, `End_position`, `Variant_Classification`, `Variant_Type`, `Reference_Allele`, `Tumor_Seq_Allele1`, `Tumor_Seq_Allele2`, `Tumor_Sample_Barcode`

Some tools in MuSiC require additional non-standard columns in the MAF (ordering doesn't matter, as long as they're after column 34). For example, the tools `cosmic` and `proximity` need the MAF to specify the `transcript_name` (doesn't matter what convention as long as its consistent), `amino_acid_change`([protein change in HGVS][3]), and `c_position`([codon change in HGVS][4]) that each variant was annotated to.

[1]: https://wiki.nci.nih.gov/pages/viewpage.action?pageId=26187384
[2]: https://github.com/genome/gms-core/blob/master/lib/perl/Genome/Model/Tools/Music/Bmr/CalcBmr.pm#L363
[3]: http://www.hgvs.org/mutnomen/recs-prot.html
[4]: http://www.hgvs.org/mutnomen/recs-DNA.html
Question: Of course, as a bioinformatician, I am aware of many large-scale open-source bioinformatics datasets, such as

- The ENCODE consortium (www.encodeproject.org, RNA-Seq, ChIP-Seq and so on),
- The Roadmap Epigenomics consortium (www.roadmapepigenomics.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The IHEC consortium (www.ihec-epigenomes.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The TCGA/ICGC consortia (www.cancergenome.nih.gov, www.icgc.org, large-scale cancer data, DNA-Seq, RNA-Seq, etc.) and
- The LINCS consortium (www.lincscloud.org/l1000 , gene expression for more than a million of different perurbation experiments).

I am wondering, however, what other wonderful datasets, the are both large and open-source, are currently available. That might include things like RNA-Seq, Chip-Seq, Bisulfite-Seq, whole genome sequencing, WGAS, and many other assays (not necessarily NGS-related, though that is what I am mostly looking for).

Also things like the (neural) connectome of certain species (in any event large data) could be of interest.

There are quite some GEO datasets that at least partially fulfill these requirements, but most are simply having to few data samples in order to be interesting to me.

Your comments are greatly appreciated!

Answer: Loss of function phenotypes: [GenomeRNAi][1].

[1]: http://www.genomernai.org/
Question: Hi all

I'd like to know if anyone has or might know where to get a list of SNPs in high LD for the last human genome build (GRCh37). I'm currently working with the Brazilian population. I'm trying to do PCA for population structure analysis and, therefore, need to remove high LD regions from my dataset to get a low biased PC axes.

Thanks a lot in advance for the attention!

Best regards,
Pedro.
PS: Just --indep-pairwise option in PLINK won't solve my problem

Answer: The LD data for hapmap / hg18 are available here: http://hapmap.ncbi.nlm.nih.gov/downloads/ld_data/?N=D

you could run liftOver to map those positions to hg19.
Question: Dear all,

I am trying to find a R package to plot LD blocks.

I found LDheatmap package but I would prefer a graphical representation as in haploview with LD blocks, and options to add lines or select markers on the graph.

Thank you for your help.

Steven

Answer: If you are particular about an R package, I would like to suggest snp.plotter package. See snp.plotter reference manual here.

Example plot from snp.plotter:

In addition, I would like to suggest LocusZoom or SNAP for visualization of LD blocks.

There is a web-interface and standalone package available for LocusZoom.

Example figure using LocusZoom:

Example figure using SNAP:

Manuscript describing snp.plotter here, LocusZoom here and SNAP here.
Question: Hi all,

I have been using R to create and manipulate phylogenetic trees for a while, but I was wondering what the differences are between the rooting functions in different packages, `ape::root`, `ggtree::reroot` and `phytools::reroot`.

I have attached an example. Every tree looks slightly different, and the tree created by `ggtree::reroot` seems best to me. Could anyone explain the differences?

tree <- read.tree(text = "(t6:0.3806255851,(t7:0.5440872659,(t5:0.6203179485,t1:0.7089423968):0.3421749519):0.2618428892,((t4:0.5886056426,t3:0.8614832051):0.09094934678,(t8:0.6766210056,t2:0.1412485428):0.4451906278):0.2356684168);")

par(mfcol=c(2,2), mar=c(0.5,0.5,1.3,0.5))

plot(tree)
nodelabels()
title("original")

tree1 <- ape::root(tree, node = 11)
plot(tree1)
nodelabels()
title("ape::root")

tree2 <- ggtree::reroot(tree, 11)
plot(tree2)
nodelabels()
title("ggtree::reroot")

tree3 <- phytools::reroot(tree, 11)
plot(tree3)
nodelabels()
title("phytools::reroot")

![Example][1]

Any explanation would be greatly appreciated.

Best wishes,

Michael

[1]: https://i.imgur.com/WM5XPgf.png

Answer: ggtree::reroot is *not* re-rooting the tree at node 11, but on the edge connecting nodes 10 & 11, so is incorrect since the position used along this edge (its midpoint?) is arbitrary. If desired this can be done using phytools::reroot by specifying the position at which the tree should be re-rooted. (I suspect it can also now be done by ape::root.)
Question: I wonder if an annotation database exists which contains all known feature across the whole genom hg19 using keywords from [sequence ontology][1] ?

[1]: http://www.sequenceontology.org/

Answer: You can use Variant Effect Predictor or SnpEff to annotate Sequence Ontology features.

http://www.ensembl.org/info/genome/variation/predicted_data.html

http://snpeff.sourceforge.net/
Question: I like to have exactly 1 match reads (6729690) from output of bowtie2 aligner. the blow cod is what I have used for aligning.

bowtie2 --very-sensitive --score-min C,0,0 -p 8 -x Genome_Index input.fastq.gz -S output.sam

8506903 reads; of these:
8506903 (100.00%) were unpaired; of these:
923759 (10.86%) aligned 0 times
6729690 (79.11%) aligned exactly 1 time
853454 (10.03%) aligned >1 times
89.14% overall alignment rate

Also, I used :

samtools view -hf "AS:i:" output.sam | grep -v "XS:i:" > unique_output.sam
samtools view unique_output.sam | wc -l

but I got 7653449 reads that is (923759 + 6729690)

Is my code right?

How can I get only 6729690 reads?

Answer: You want to remove unaligned reads ?

samtools view -F 4
Question: I've got a few thousand small bam files produced against the exact same reference, and I want to merge them into one single big bam file. What is the best way to do that?

Should I do this iteratively or can I pass a long list of bam files to samtools/picard/etc in one go?

Edited, since this is now partially solved. In my terminal, the methods below works for up to 4092 files. More than that raises an error:

samtools merge all.bam *.bam
samtools merge all.bam `find /basedir/ -name "*myfiles*.bam"`
samtools merge all.bam /basedir/*/???/*myfiles*.bam

Answer: Both are possible, but the second is probably best and easier

Put all your bam files in a folder (or create a folder with symbolic links to all sam files you want merge) and then

samtools merge finalBamFile.bam *.bam

See [here][1] for options about samtools (header and so on)

Editing to address the issue with many file.

You can take a three steps approach.

First you print the header

samtools view -H bar/foo/anyFile.bam > allSeq.sam

then you append sequences from all files using find

find ./ -name \*.bam -exec samtools view {} \; >> allSeq.sam

and finally convert to bam

samtools view -bh allSeq.sam > allSeq.bam

I can't quite test it at the moment, but should work

[1]: http://samtools.sourceforge.net/samtools.shtml
Question: I was comparing the percent of FLT3 genes in AML from TCGA data with the mutations found in our exome data. Our exome data only reveals less than 50% of FLT3 mutations compared to TCGA samples.
Does anyone know what capture method TCGA consortium used for AML samples? Is there any reason why we can't see half of FLTs in our cohort compared to TCGA? Also, what could be the reasons that all mutations in FLT3 can't be determined by exome methods?

Answer: tl;dr: You should definitely not expect to recover all FLT3 mutations from exome sequencing and stock variant calling.

I was the lead analyst for the genomics in the AML TCGA paper. Finding all of those FLT3 ITDs was a giant pain in the ass. At the very least, you should also use the targeted validation data. Even, then, I wouldn't expect to find everything. We did some targeted 454, and even some Sanger, I think. Even if you have good coverage of the region, variant callers generally suck at finding the ITD, because of the short reads and repetitive sequence. Many required manual inspection to resolve.

There are some additional tips in this previous thread: https://www.biostars.org/p/158375/ Specifically, I'd have someone sit down with every sample and pull up that region, looking carefully for soft-clipping, etc that will help you identify the event. Good luck!