cannin/biostars_qa huggingface dataset:
Question: Hello Everyone,

We have recently generated two de-novo transcriptomics assembly for two different but related species. These new transcripts seem quite good on the basis of quality measurements, completeness and alignment with the previously sequenced genome and annotation. We were able to pick up novel genes and previously unannotated transcripts. In order to pick the alternate spliced transcripts we are currently trying running PASA (Program to Assemble Spliced Alignment) annotation.

After PASA annotation, when I made a comparison between trinity transcripts, already existing annotation and PASA annotation in IGB. I find out a case where trinity transcripts fully supported by previously curated annotation as well as the RNASeq data, but no PASA annotation. PASA annotation only shows one fragment, rest of the parts not even present in the valid and failed .gff3 file generated by PASA [Figure]. This leads to few questions - for which I don't have any answer. Here are my questions:

1. Why there is this difference in the PASA annotation- when there is already manual curation and RNAseq depth evidences are present? And which one is correct?
2. PASA annotation using blat and gmap for the alignment of transcripts to the genome. And we have also used blat to align the transcript to the genome. Then why there is different in two blat alignment?

**Data used for Transcriptomics assembly** = 100 bp Paired end reads, non-strand specific

**Attached figure description**:-

- Dark Blue colour = New assembled transcript annotation [this is the annotation generated by aligning assembled transcripts with the genome using blat].
- Orange = existing curated annotation.
- Red = PASA annotation
- Blue colours = shows valid alignment annotation for blat and gmap respectively.
- Read Depth = Green for control and Dark red = Knock-out

Figure: ![< image not found >][1]

I am very much looking forward for the reply. Any suggestions/view would be very helpful.

Many Thanks
Reema,

PS = I also posting the same post in the seq-answer.

[1]: http://www.compbio.dundee.ac.uk/user/rsingh/Figure.jpeg

Answer: Hi Reema,

Our senior software developer was able to look over your .gff3 file. The issue has to do with how IGB currently displays .gff3 data. It's something we should be able to fix with a small patch.

The .bed file does display correctly, so if you have a choice in the output, go with the .bed file for now.

Thank you for helping us identify this issue!

Nowlan
Question: Hi all,

I am annotating my new plant genome now and am working with Maker and its very detailed tutorial (http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/MAKER_Tutorial_for_WGS_Assembly_and_Annotation_Winter_School_2018) . I have read a few really helpful posts about Maker here as well, but i still have some questions.

1. SNAP training. How do you actually know that it is enough to train and you can run your final Maker run? I have tried to run it several time and there is a difference in the number of genes every time. It is actually a kind of sinusoidal graph - number of genes are going up and down... So when do you stop? Or how do you know that SNAP is trained? Do you wait until the plateau? How many times did you do the training and why?

2. My genome has unusually high repeat content. This is why I decided to create its own repeat library with repeatModeler. The question is where on the option file do I add this repeat library?

###THANKS a lot for your help,
Alex

Answer: You can specify a custom repeat library (in FASTA format) with `rmlib` in the Repeat Masking section of the `make_opts.ctl` file

#-----Repeat Masking (leave values blank to skip repeat masking)
model_org=all #select a model organism for RepBase masking in RepeatMasker
rmlib=repeatlibrary.fa #provide an organism specific repeat library in fasta format for RepeatMasker
repeat_protein=/opt/maker/data/te_proteins.fasta #provide a fasta file of transposable element proteins
rm_gff= #pre-identified repeat elements from an external GFF3 file
prok_rm=0 #forces MAKER to repeatmask prokaryotes (no reason to change this), 1 = yes, 0 = no
softmask=1 #use soft-masking rather than hard-masking in BLAST (i.e. seg and dust filtering)`

This is an example Repeat Masking section

You might also consider running ProtExcluder on the output of RepeatModeler

http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/Repeat_Library_Construction-Basic

# Run blastx then ProtExcluder to excluce known protein sequences from RepeatModeler library
/usr/bin/blastx -num_threads 75 -db /genetics/elbers/maker/uniprot_sprot.fasta -evalue 1e-6 \
-query repeatlibrary.fa -out repeatlibrary.fa.blast

/opt/ProtExcluder1.1/ProtExcluder.pl -f 50 repeatlibrary.fa.blast repeatlibrary.fa
# output of ProtExcluder is "temp"
# rename temp to whatever you desire
mv temp repeatlibrary.fa2
Question: Hi,

I used data set from Encode consortium for my package development, due to size of actual peak files are rather big, I can't use these data set for my package use. Because actual size of package resulted from R CMD build must be less than 4Mb on disk, I have to use rather small peak file as an example data for my package . In Encode sample's data set, each peak files contains around 100,000 peaks each. How can I edit rather big bed files in order to keep particular chromosome ? Is there any handy tools to edit peak files ? Thanks in advance :)

Best regards :

Jurat

Answer: If you have GNU Parallel installed, you can use this with [BEDOPS `bedextract`][1] to very quickly split a BED file by chromosome:

$ bedextract --list-chr input.bed | parallel "bedextract {} input.bed > input.{}.bed"

You can then use my [`sample` utility][2] or GNU `shuf` to uniformly sample without replacement:

$ sample -k ${SAMPLE_SIZE} input.chrN.bed > input.chrN.sample.bed

Or:

$ shuf --head-count=${SAMPLE_SIZE} input.chrN.bed > input.chrN.sample.bed

[1]: http://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedextract.html
[2]: https://github.com/alexpreynolds/sample
Question: Currently when calling variants you have to call variants on all people simultaneously, which allows for rare variant positions being check and called in every person without such variant.

My question is, why wouldn't it make more sense to just determine variant in every genomic position (3.2Billion) and that's it? You could merge this dataset with others easily, moreover, no need to do joint variant calling etc. It is very frustrating to see such an inefficiency becoming default.

Answer: You are correct, but you have conflated two issues - variant calling and warehousing. First, I don't think you will see joint genotyping being routinely done in 5 years - single sample calling with instrument-specific training sets are where things are headed. Secondly, the bigger groups are moving or have already moved to variant warehousing rather than VCFs - TileDB, Google Variant Transforms, GenomicsDB - and these want gVCFs as input.
Question: Hi,
Suppose I have a bam file and a vcf file containing variant calling result. I want to extract only reads with their mate that support variant allele in the vcf. It would be nice to get those reads in bam format. I tried googling such tools to do this and found like VariantBAM but it reports both reads that supporting and not supporting variant.

Thanks

Answer: I quickly wrote something : http://lindenb.github.io/jvarkit/Biostar322664.html

caveat: BAM files must be sorted with picard **SortSam/queryname** and variants are loaded in memory and only SNP are considered.

$ java -jar picard.jar SortSam I=src/test/resources/S1.bam O=query.bam SO=queryname
$ java -jar dist/biostar322664.jar -V src/test/resources/S1.vcf.gz query.bam

(...)
RF02_358_926_2:0:0_2:1:0_83	83	RF02	857	60	70M	=	358	-569	GACGTGAACTATATAATTAAAATGGACAGAAATCTGCCATCAACAGCTAGATATATAAGACCTAATTTAC	2222222222222222222222222222222222222222222222222222222222222222222222	RG:Z:S1	NM:i:3	AS:i:55	XS:i:0
RF02_362_917_2:0:0_2:1:0_6f	147	RF02	848	60	70M	=	362	-556	ATAAGGAATCACGTTAACTATATACTTAAAATGGACTGAAATCTGCCATCAACAGCTAGATATATAAGAC	2222222222222222222222222222222222222222222222222222222222222222222222	RG:Z:S1	NM:i:3	AS:i:55	XS:i:0
(...)

it's not fully tested, tell me if something is looks wrong.
Question: Hi,

I have a barplot like this:

![enter image description here][1]

[1]: https://imgur.com/pizrVlq.jpg

I created it with this code:

toplot.N <- data.frame( set=c("FLCN", "All SNPs", "eQTL from 103 genes","All eQTL"),
FDR =c(FDR1FI,FDR2FI,FDR3FI,FDR4FI))

p<-ggplot(toplot.N, aes(x=set, y=FDR, fill=set)) + labs(y = "pi_1")+
geom_bar(stat="identity")+theme(axis.title.x = element_blank())+
geom_hline(yintercept=0.05, linetype="dashed", color = "red")
p

How would I change this so that the order of bars is: "All SNPs","All eQTL","eQTL from 103 genes","FLCN"

Thanks
Ana

Answer: We could change order of factor levels before plotting:

toplot.N$set <- factor(toplot.N$set,
levels = c("All SNPs", "All eQTL", "eQTL from 103 genes", "FLCN"))
Question: hi
can anyone tell me the name of the software for performing the alignment and constructing the phylogenetic tree of whole genome. thanks in advance.

Answer: I am not aware of an easy way to construct reliable species trees based on complete genomes. The general approach that you need to take is to pick one or more genes based on which to base your phylogeny. This could be either 16S rRNA, all ribosomal-protein-coding genes, or other highly conserved genes that are universally present and rarely subject to gene duplications or lateral gene transfer.

Once you have picked the genes, you need to make a multiple sequence alignment(s). You need to do this for each of the genes that you want to use for your phylogeny. For this I would tend to use either [muscle][1] or [mafft][2]. After that I would use [Gblocks][3] to extract the conserved blocks in the alignment(s) in order to not use potentially misaligned parts as the basis for tree building.

If you decided to use multiple genes as the basis for your phylogeny, you now have to make a big decision, namely whether to go for a concatenated alignment approach or a supertree approach. In the first case, you would concatenate all of the multiple alignments and use the resulting big alignment as input for a phylogenetic tree reconstruction program, for example [PhyML][4]. In the second case, you would use such a program to make a separate tree for each of the genes of interest, and subsequently use one of several supertree programs to derive a consensus tree based on these. If you went for just using a single gene as the basis for your tree, you obviously just build a tree for that one gene and you are done.

I hope this helps, although it is certainly very far from a "push of a button" solution.

[1]: http://www.drive5.com/muscle/
[2]: http://mafft.cbrc.jp/alignment/software/
[3]: http://molevol.cmima.csic.es/castresana/Gblocks.html
[4]: http://atgc.lirmm.fr/phyml/
Question: Why is it so difficult to make things in ggplot2 , i like the way it helps in customisation but the curve is steep nevertheless

Here is my sample dataframe

df <- gene HSC CMP
ENSG00000158292.6 1.8102636 2.456869
ENSG00000162496.6 2.6796705 6.203838
ENSG00000117115.10 3.4509115 5.555739
ENSG00000159423.14 3.6809277 5.063446
ENSG00000053372.4 5.7089974 6.851090

If i have plot a boxplot i can simply write this `boxplot(df[,-1],col=c("red","blue"))`

I get a boxplot but when im trying with ggplot2 im having difficult time

ex <- melt(df, id.vars=c("HSC", "CMP"))
ggplot(data = ex,
aes(x = CMP, y = HSC)) +
geom_boxplot()

I get a single boxplot what i want is i get a box plot for HSC and CMP which i got when i use simple base R boxplot .

Any help or suggestion would be highly appreciated with my ggplot2 code

Answer: Devon got there before me but as he mentioned the id.vars needs to be set to 'gene'

Here's a boxplot with scatterplot overlay for anyone else arriving here from Google.

I do agree that ggplot can be difficult to work with. Many functions redundant in the sense that they do the same thing as other but have different names, and conflicts frequently arise. That said, if you can master ggplot, then you can produce very nice graphics for publications.

require(reshape2)
require(ggplot2)

ex <- melt(df, id.vars=c("gene"))
colnames(ex) <- c("gene","group","exprs")

ggplot(data=ex, aes(x=group, y=exprs)) +

geom_boxplot(position=position_dodge(width=0.5), outlier.shape=17, outlier.colour="red", outlier.size=0.1, aes(fill=group)) +

#Choose which colours to use; otherwise, ggplot2 choose automatically
#scale_color_manual(values=c("red3", "white", "blue")) + #for scatter plot dots
scale_fill_manual(values=c("red", "royalblue")) + #for boxplot

#Add the scatter points (treats outliers same as 'inliers')
geom_jitter(position=position_jitter(width=0.3), size=3.0, colour="black") +

#Set the size of the plotting window
theme_bw(base_size=24) +

#Modify various aspects of the plot text and legend
theme(
legend.position="none",
legend.background=element_rect(),
plot.title=element_text(angle=0, size=14, face="bold", vjust=1),

axis.text.x=element_text(angle=45, size=14, face="bold", hjust=1.10),
axis.text.y=element_text(angle=0, size=14, face="bold", vjust=0.5),
axis.title=element_text(size=14, face="bold"),

#Legend
legend.key=element_blank(),	#removes the border
legend.key.size=unit(1, "cm"),	#Sets overall area/size of the legend
legend.text=element_text(size=12),	#Text size
title=element_text(size=12)) +	#Title text size

#Change the size of the icons/symbols in the legend
guides(colour=guide_legend(override.aes=list(size=2.5))) +

#Set x- and y-axes labels
xlab("Stem cell class") +
ylab("Expression") +

#ylim(0, 0) +

ggtitle("My plot")
Question: Hello!

I currently have a 255k assembled contigs transcriptome (de novo, Trinity) and I want to retrieve indentifiers and sequences in fasta format of the contigs that had hits in BLASTx (around 50%). There's any way to filter using the BLASTx output, an script and the transcriptome itself? O I should search an approach of parsing BLASTx output to contain sequences?.

Thanks in advance!

Answer: Many people have encountered the problem before. See their ports:

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

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

or this guide:

http://sfg.stanford.edu/BLAST.html
Question: I have 4 control and 4 treated samples of RNA seq data, for generating gene expression network should i consider adding the normalised counts of control samples in the gene expression matrix?

Answer: Yes, you would tend to include all of the group when making a coexpression network. You need larger sample numbers to make a reliable network anyway, so you can only achieve that by including all of your samples.
Question: I am looking at RNA-seq data, which I have little experience in. I notice that for many genes, there are reliable alignments (i.e. with high mapping quality) to introns. I understand that some of them are due to unannotated transcripts, but in many regions, this does not seem to be the major cause. The intronic read hits do not seem to be purely caused by alignments artifacts, either, because the pattern is tissue specific (though this is not a compelling evidence). Another possible explanation is that this observation is due to noisy transcripts (Pickrell et al, 2010), but this seems to be a big effect: for some long genes, there are far more intronic hits than exonic hits.

I guess those who study RNA-seq data must have noticed the intronic hits for years. What is cause of the large amount of intronic read hits? Is it caused by alignment/library prep artifacts or noisy transcription? Are there papers addressing this? Thanks.

EDIT: my conclusion. I was looking at ERR030882 from Illumina BodyMap (brain). The sample were processed with oligo-dT. I am using the gencode exon annotations, including all the pseudogenes, lincRNA and known processed transcripts, totalling ~112Mbp. The initial analysis reveals ~80% of bases mapped to exons. Nonetheless, if I only look at read pairs with insert size larger than 311bp (~10% of the original data), 98.2% of these spliced read pairs are mapped to known exons, suggesting that the vast majority of the intronic and intergenic read pairs are unspliced. It is possible that some unspliced pairs come from unknown single-exon transcripts with intact polyA tail, but contaminations seem the leading cause overall.

Answer: Two simple reasons:

1) Genomic DNA has contaminated the RNA-Seq sample, likely at the mRNA isolation step. This would look like sequence data from both strands of the intron.

2) There is unspliced mRNA in the sample. This would give data for the strand that encodes the gene, but within that intron.

There could be other explanations, but these are two principle ones that come to mind.
Question: Hi,

For the standalone version, can you answer the following please:

1. How does OMA deal with jobs that do not finish in one submission i.e. when working on clusters jobs can only run for a certain amount of time. To 'continue' the job do I simply resubmit the job and OMA detects where it was from the Cache?

Thank you,
R

Answer: In the Cache/AllAll/ directory, all files that are gzipped represent job chunks that have successfully completed and will be used. The files that are not gzipped represent job chunks that were being processed when your job died. You can delete these and they will be restarted next time you run OMA.
Question: Hi all,

In the article of Miwa: https://www.ncbi.nlm.nih.gov/pubmed/15816807
They mention a polymorphism M429V of the gene ABCG8, I need the rs number of this SNP but I couldn't find it. Could someone provide this SNP.

Thank you in advance.

Answer: Using [backlocate][1], it looks like the variant would be in the BED interval 2:44100998-44101000 .

#User.Gene	AA1	petide.pos.1	AA2	transcript.name	transcript.id	transcript.strand	transcript.AA	index0.in.rna	wild.codon	potential.var.codons	base.in.rna	chromosome	index0.in.genomic	exon	messages	extra.user.data
ABCG8	Met	429	Val	ABCG8	ENST00000272286	+	M	1284	ATG	GTG	A	2	44100998	ENST00000272286.Exon9	.	.
ABCG8	Met	429	Val	ABCG8	ENST00000272286	+	M	1285	ATG	GTG	T	2	44100999	ENST00000272286.Exon9	.	.
ABCG8	Met	429	Val	ABCG8	ENST00000272286	+	M	1286	ATG	GTG	G	2	44101000	ENST00000272286.Exon9

.	.

This position would be https://www.ncbi.nlm.nih.gov/snp/rs147194762

[1]: http://lindenb.github.io/jvarkit/BackLocate.html
Question: Hi everone,

I searched Biostars for BAM/SAM to FASTA conversion method, and found the tools EMBOSS Picard could do this (https://www.biostars.org/p/6970/). I am wondering is there any perl script to accomplish this work? I should note that my seq data is strand-specific seq, so the command `samtools view filename.bam | awk '{OFS="\t"; print ">"$1"\n"$10}'` may not be very good.

I appreciate any of your solutions. Thanks a lot!

Answer: For future record: samtools versions >1.3 can convert bam to fasta directly via `samtools fasta`

samtools --help

Program: samtools (Tools for alignments in the SAM format)
Version: 1.4 (using htslib 1.4)

Usage: samtools [options]

Commands:
[...]
-- File operations
[...]
fastq converts a BAM to a FASTQ
fasta converts a BAM to a FASTA
[...]
Question: Dear all,

please could you advise, what is a typical range of AF (allele fractions) of somatic mutations in the cancer samples in order to consider them clonal or sub-clonal ? Thank you ...

-- bogdan

Answer: In addition to Kevin's answer, [Quantification of subclonal selection in cancer from bulk sequencing data][1] and references therein may be another useful read.

If you want some sensible cutoff for some rough analysis or filtering, I would say 20% AF is about right, but again, as Kevin says it depends...

[1]: https://www.ncbi.nlm.nih.gov/pubmed/29808029
Question: I am trying to interpret the component graphs from my Trinity run. I rendered a couple of graphs of the components (using c*.graph.out files) from a Trinity assembly and noticed that some components had a structure where the root node (-1 node) is in the middle of a "linear" sequence of nodes.

I uploaded my ipython notebook here with 3 of the graphs I rendered:
http://nbviewer.ipython.org/github/damiankao/trinity-visualization/blob/master/trinity_vis.ipynb

The first component (c445) looks normal to me, with a root node (in red) that connects to one linear sequence and eventually splitting into two branches followed by a merge that could possibly indicate isoforms.

But the second and third component graphs showed the root node in the middle of a "linear" region. Furthermore, for the second component graph, the probable paths were the two "arms" of the root node.

There are no shared k-1-mers between the nodes on either side of the root node in the second and third component. How does the bundling of contigs work in this case? Is it putting these contigs together based on pair-end reads? And what exactly is the -1 root node?

Answer: I just got a response from the Trinity mailing list via Brian Hass:

> The -1 is the root node for the de Bruijn graph. A way you can end up with multiple 'arms' in the graph like that is if there are multiple inchworm contigs that are clustered together based on paired-read links (from the bowtie alignment step). This way, they end up having the same 'component' number in the accession string (ie. (`c\d+`) of the `c\d+_g\d+_i\d+` accession naming format. This often happens when transcripts for a given gene are fragmented.
Question: I have Total RNA TrueSeq Illumina Stranded library (human). My goal is to find novel (and non-novel) non-coding transcripts in my data (experimental vs control).

After a LOT of Google-fu and asking questions on this website, this is the methodology that I am currently using -

1. Align the fasta files with STAR to hg38
2. Assemble transcripts for each sample, merge transcripts from all samples (to get a unified transcriptome that represents all the samples), and estimate transcript abundances - all using Stringtie (protocol paper - https://www.nature.com/articles/nprot.2016.095#procedure)
3. Use tximport to infer integer counts from the Stringtie transcript abundances and export it to DESeq2.

I wanted to know if this methodology makes sense. Is there anything for which a better method makes more sense.
I hope my question is not too broad, given that I do specify the exact pipeline I am employing :)

Answer: Hi, you can follow the below given pipeline for complete analyses

1. Algin reads to the genome using the STAR (Already completed)
2. Run stringtie to assemble transcripts with default parameters which (filter out transcript below 1 FPKM)
3. Merge all transcriptome assemlby into 1 merged assembly ( stringtie --merge)
4. Extract protein coding genes in gtf format from ENSEMBLE GTF
5. Download already known noncoding RNAs ( ENSEMBLE + LNCIPEDIA)
6. Merge all known noncoding annotation in to 1 gtf file by ( using cuffmerge)
7. gffcompare the protein coding and all knonw protein codng gtf to your assembled merged gtf file.
8 Extract the classes "u", "i" and "X" with bash script (For more details check gffcompare classes.)

9. Now extract the MSTRG id and extract gtf file for this using the grep or something.
10 Extract fasta file for novel noncoding using gffread
11. Filter the transcript wiht length less then 200 and single exons transcript.
12. Use CPAT and PLEK to predit the coding potenial for the novel transcripts using the transcript sequence.
13. Further filter the the transcritp by using NCBI blast with known protein homology based ( use only 3 frames)
14. This would be your final list for the transcripts.
15 Rerun the stringtie with -e -B option to get count in FPKM, Coverage and TPM.

16 Run DESEQ2 For differential expression

17. Annotate the novel transcript by looking into near by genes by bedoops closet tool.
18. Further find out signifigance by creating the corrleation beteween protein coding and your novel transcripts.

19 Gene enrchimnent anlaysis of the near by genes and protein coding gens for understaing mechanism