cannin/biostars_qa huggingface dataset:
Question: The plot below is what I want to plot:
![enter image description here][1]

Basically, the plot contains the log(Pi1/Pi2) values for specific positions on each chromosome. (It looks like a Manhattan plot but with negative values)

I have the following data frame

Chromosome Position Log_Pi

1 102331 -0.552
2 220231 2.32
5 233433 -2.200
.
.
.
7 522631 -1.023
7 322512 1.2231
10 233356 -0.223
12 5666932 0.2235

I have tried:

ggplot(data = Pi_PLOT, aes(x = Chromosome, y = log_pi)) + geom_point(data = Pi_PLOT, aes(x = Position, y = log_pi, fill = Chromosome)) + labs(x = "Chromosome", y = "log_pi")

[1]: /media/images/b3c1418c-55f8-4613-9257-e700c256

Answer: I figured out by using the [CMplot package](https://cran.r-project.org/package=CMplot) in R.

Thanks
Question: Hi,

1- In the table ´Normalization method' [here][1] says that CPM (counts per million) can be used for gene count comparisons between replicates of the same sample group.

1.1 Does it mean that for eg. we can compare one gene from a sample of group 'control' with same gene of another 'control' sample but that we cannot compare a gene in 'control' sample with same gene in a 'treatment' sample?

1.2 If so, then when looking for a heatmap with CPM values cannot we for e.g. identify genes that seem to have a higher expression in 'treatment' samples than in 'control' samples? Do we need to use a different normalization method?

2- In same table says that CPM cannot be used for within sample comparisons.

2.1 Does it mean we cannot compare different genes of the same sample?

2.2 What if when looking to CPM heatmap it seems one gene is varying more between 'control' and 'treatment' than the other. Can we make this conclusion if heatmap plots CPM values?

[1]: https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html

Answer: As recommended in this presentation, I would not use per-million methods for anything as there are better methods now. Check [this][1] video to get an idea why per-million based methods are not optimal and [this][2] one on how the normalization in e.g. DESeq2 works.

Towards your questions:

1 - you can use it but it is not recommended for DE analysis, so better don't use it at all

1.1 - Simply normalize the entire dataset with edgeR or DESeq2 and do comparisons with these values

1.2 - do not use CPM values for a heatmap, use logged/normalized counts, like those produced by the `vst` or `rlog` functions in DESeq2. Using non-log counts will bias the heatmap towards highly expressed genes. These video series I inked above also have a video about logs in case you care.

2 - true, because it does not normalize for gene length, so longer genes inherently have higher counts than short genes.

2.1 - one probably could, but not without adjusting for gene length (use the search function on this, there are plenty of posts on that matter already out there).

2.2 - it might give you an idea but you should use appropriate statistics to infer differentially expressed genes.

[1]: https://statquest.org/2015/07/09/rpkm-fpkm-and-tpm-clearly-explained/
[2]: https://www.youtube.com/watch?v=UFB993xufUU
Question: **Can I transfer this first format to the second one just by basic shell procession or awk or sed on linux? This is a toy example:**

This kind of text file is what I have, three cols, col2 and col3 like range, left close and right open,
chr1 0 2 0
chr1 2 6 1.5
chr2 0 3 0
chr2 3 10 2.1

Transfer to describe each position as:
chr1 0 0
chr1 1 0
chr1 2 1.5
chr1 3 1.5
chr1 4 1.5
chr1 5 1.5
chr2 0 0
chr2 1 0
chr2 2 0
chr2 3 2.1
...
chr2 9 2.1

Someone has idea how to solve? Thanks!!!

Answer: Try:

awk -F "\t" -v OFS="\t" '{while ($2 < $3) { print $1,$2,$4; ++$2 }}' test.txt

or

awk -F "\t" -v OFS="\t" '{for ($2=$2; $2<$3;$2++) {print $1,$2,$4}}' test.txt

output:

chr1	0	0
chr1	1	0
chr1	2	1.5
chr1	3	1.5
chr1	4	1.5
chr1	5	1.5
chr2	0	0
chr2	1	0
chr2	2	0
chr2	3	2.1
chr2	4	2.1
chr2	5	2.1
chr2	6	2.1
chr2	7	2.1
chr2	8	2.1
chr2	9	2.1
Question: Hi,

I was using the pheatmap to cluster all the samples and do the visualization, but i was confused that how to increase the size legend which shows like that.

Thank you so much.

Answer: Take a back up of your code posted above, and try plotting with this basic code (OP code). If you can see the legend, then keep optimizing heatmap. I think unoptimized cell and other parameters are issue here in code (as data is not furnished) :

pheatmap(data,
method = c("pearson"),
clustering_method = "complete",
treeheight_row = 40,
treeheight_col = 40,
cluster_row = T,
cluster_col = T,
show_rownames = T,
show_colnames = T,
legend = T,
fontsize = 10,
color = colorRampPalette(brewer.pal(9,"Reds"))(400))
Question: Hi,

I have several bed files for which I would like to get conservation information. I have downloaded the 100wayPhastCons.bw file from UCSC and my plan is to now get the conservation scores for my bed intervals from this file.

Is there a way to achieve this ?

Thanks,

Answer: Here's a pipeline you can use:

1. Convert the conservation bigWig (.bw) to Wiggle (.wig) with `bigWigToWig`, using `/dev/fd/1` as a filename for standard output. This should allow you to pipe the result to BEDOPS `wig2bed` and get a sorted five-column BED file:

$ bigWigToWig conservation.bw /dev/fd/1 | wig2bed - > conservation.bed5

2. Split your intervals up into per-base elements with `bedops --chop 1`.

$ bedops --chop 1 intervals.bed > perbase.bed3

3. Run `bedmap` to map your per-base elements to the BED-formatted conservation signal. Normally, unmapped elements report no value. You can use the following command to replace unmapped bases with a score of `0`:

$ bedmap --count --echo --echo-map-score --delim '\t' perbase.bed3 conservation.bed5 | awk '{ if ($1=="0") { print $0"0.0"; } else { print $0; } }' | cut -f2- > answer.bed

This will only make sense if a missing conservation score should be zero. You could alternatively replace `0` with `NaN` or a similar IEEE-ready placeholder for a non-number.

Finally, if you want to put everything into a one-liner, use file substitutions:

$ bedmap --count --echo --echo-map-score --delim '\t' <(bedops --chop 1 intervals.bed) <(bigWigToWig conservation.bw /dev/fd/1 | wig2bed - ) | awk '{ if ($1=="0") { print $0"0.0"; } else { print $0; } }' | cut -f2- > answer.bed

If you're going to do repeated queries with different sets of intervals, it probably makes sense to do the bigWig-to-BED5 conversion step only once, and drop that in as `conservation.bed5` where you need it in mapping.
Question: Dear all,

I am trying to convert a genbank file into genome feature format version 3. Using the [GFF3 online validator][1] the file I have prepared are OK but there are issues with some of the phases of the CDS. I am summarizing the issues in this table:

entry	phase	start	end	delta	Delta third
1	0	76014	76437	423	141.00
2	2	76532	76794	262	87.33
3	0	76901	77296	395	131.67
4	0	82169	82217	48	16.00
5	2	82326	83644	1,318	439.33
6	0	83960	85309	1,349	449.67
7	1	86174	88442	2,268	756.00
8	0	88544	89010	466	155.33
9	1	89700	90945	1,245	415.00
10	0	91042	91496	454	151.33

entry	phase	Minus 1	Minus 1 third	Minus 2	Minus 2 third
1	0	422	140.67	421	140.33
2	2	261	87.00 260	86.67
3	0	394	131.33	393	131.00
4	0	47	15.67	46	15.33
5	2	1317	439.00 1316	438.67
6	0	1348	449.33	1347	449.00
8	1	2267	755.67	2266	755.33
9	0	465	155.00 464	154.67
10	1	1244	414.67	1243	414.33
11	0	453	151.00 452	150.67
`entry` is a given CDS feature.

`phase` is the suggested phase from the validator

`start` is the start position of the CDS and `end` its end position

`delta` is the difference `end - start`

`delta third` is `delta/3`

`Minus 1` is `delta -1` and `Minus 1 third` is `(delta -1)/3`

`Minus 2` is `delta -2` and `Minus 2 third` is `(delta -2)/3`

Taking the entry #1, its length delta is divisible by 3 so it makes sense that the validator accepted a phase of 0. Same thing for feat. 4, 6, 8, 10.

The second entry it is not divisible by 3 directly so it is understandable that the validator has flagged it out. Shortening the length of the feature by 1 nucleotide (Minus 1), the feature is now divisible by 3, thus I expected to change the phase to 1, not 2. Same thing for feat. 5.

The third feature is also not divisible by 3, yet the validator did not flag it. It is divisible by 3 after removing two nucleotides, thus I thought the phase should have been 2.

Features 7 and 9 are divisible by 3, yet are flagged with a phase of 1.

It is pretty confusing for me and I haven't found many tutorials online on the subject.

How do I calculate the phase of the CDS based on the start-end positions?

Thanks

[1]: http://genometools.org/cgi-bin/gff3validator.cgi

Answer: The actual gff3 or original genbank feature might have been useful. Typically this problem occurs only in organisms with spliced proteins, for the second and later exons of a transcript.

From your example, I guess the first three CDS entries correspond to three exons of the same gene (given the 5 kb gap). To complete the last, incomplete codon from the previous exon, you require the next codon to start *phase* number of nucleotides after the indicated position.

For this to make sense, please note that gff uses 1-based coordinates, so your feature length is `end-start+1`. Hence your first feature is 424 nucleotides long and requires a phasing of 2 for the second feature.
Question: This is bordering on off-topic, but I routinely want to quickly serve up a set of BAM files to users for viewing via IGV and the like via a simple http server. Both python and node (http-server) servers do not seem to correctly support range requests. Short of configuring and running something more heavyweight (nginx, apache, etc.), are there other simple command-line http servers that can support quickly and easily serving up a directory?

Answer: set-up a **jetty** server? see the gist below. you just have to change the path

holderPwd.setInitParameter("resourceBase", "/home/user/path/to/directory");

https://gist.github.com/lindenb/1fa88442f829d197705a
Question: Hello guys,

from your experience with chip-seq analysis what would be the most useful visual representations of chip-seq results?
I am particularly interested in your opinion on visualizations of motif discovery, annotation and differential binding analysis and your experiences with various bioinformatics tools that you' ve been using for this purpose.

Cheers,
Nemanja

Answer: If you have a lot of samples (including replicates), it can be very useful to visualize sample clustering. Examples include heatmaps with hierarchical clustering and PCA. The [DiffBind package in Bioconductor][1] provides clustering and PCA plots, as well as MA, box, and volcano plots after performing the differential analysis. You can then feed the differential sites into some of the excellent tools above.

[1]: https://bioconductor.org/packages/release/bioc/html/DiffBind.html
Question: I have this problem so i want to take out all the exon coordinates and their respective sequences .One way is I take out the coordinates of exon from a gtf file then give that as input to ucsc genome browser , i get the DNA sequences ,but it looks like i get intronic sequences too.

So any suggestion or help how to do it better way .I m not how to proceed .As always sugesiton or help would highly appreciated

Answer: I did it like that:

1- Download [refGene.txt.gz][1] and [hg19.fasta][2] from the UCSC goldenpath. ( note: convert hg19.2bit to hg19.fa using [twoBitToFa][3] )
2- Create a bed file with exon coordiniate using my [awk script][4]

// to_transcript.awk

BEGIN
{
OFS ="\t"
}
{
name=$2
name2=$13
sens = $4 =="+" ? "forward" : "reverse"
chrom=$3

split($10,exon_starts,",")
split($11,exon_ends,",")

for (i=1; i exon.bed

3 - Extract sequence using bedtools

bedtools getfasta -fi hg19.fa -bed exon.bed -fo exon.fa -name

4 - You may want to reverse sequence for gene placed on the reverse strand. Split exon.bed into and exon_minus.bed and exon_plus.bed using grep on column 5th. And do same as I described above to have exon_minus.fa and exon_plus.fa
You can then revert your sequence using :

seqtk seq -r exon_minus.fasta> exon_minus.rev.fasta

and concat both file :

cat exon_plus.fasta exon_minus.rev.fasta > all_exons.fa

[1]: http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/refGene.txt.gz
[2]: http://hgdownload.cse.ucsc.edu/goldenpath/hg19/bigZips/hg19.2bit
[3]: http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
[4]: https://github.com/dridk/codon_biais/blob/master/scripts/to_transcript.awk
Question: I have few motifs whose location on a given DNA sequence has been identified using MATCH (from TRANSFAC). MATCH provides the location of these matrices in a numerical way as motif 1 236 (+) (for an example - motif 1 is present at location 236). I want to get a graphical output by feeding the length of DNA sequence and the the position of matrices so that its easy for presentation and comparing two different sequences with almost similar motifs

Thanks in advance for any suggestions

Answer: I think this depends on the length of your sequence.
Our MSA editor can add comments to a sequence (manually or from a file), and then a "visual summary" can show 100s of kb in a single picture.
The file format takes motif name, postion, strand, colour for showing these annotations on the sequence.
You can then save this file locally (BBB is Java)
http://athena.bioc.uvic.ca/virology-ca-tools/base-by-base/

Hillary, W., S.H. Lin, and C. Upton. 2011.
Base-By-Base version 2: single nucleotide-level analysis of whole viral genome alignments, BMC Microbial Informatics and Experimentation, 1:2

William Hillary, Song-Han Lin and Chris Upton*

Brodie, R., A.J. Smith, R.L. Roper, V. Tcherepanov, and C. Upton. 2004.
Base-By-Base: Single nucleotide-level analysis of whole viral genome alignments, BMC Bioinformatics, 5, 96.

Cheers,
Chris
Question: I estimated DEGs from RNASeq using edgeR method. But the final DEG list has several genes with low counts and are having high Fold change

logFC	logCPM	F	PValue	FDR	Ctrl1	Ctrl2	Str1	Str2
Gene1	7.673	-2.658	33.203	1.27E-05	0.000118	0	0	0.738954	0.657007
Gene5	7.605	-2.951	26.996	3.36E-05	0.000255	0	0	0.674697	0.657007
Gene35	7.499	-2.072	22.724	8.47E-05	0.0034	0	0	0.642569	0.594435
Gene67	7.462	-2.151	36.849	2.49E-06	3.32E	0	0	0.642569	0.563149

Initialy, I filtered out low count genes using

keep<-filterByExpr(y,group = group)
y<-y[keep,keep.lib.size=FALSE]

But I am not sure why genes normalized count value close to zero are retained and listed as DEGs. How can I solve this? I was thinking to use min.count parameter in filterByExpr.
But, is that the right way? Suppose if I use, `keep<-filterByExpr(y,group = group, min.count=500)`, the genes with greater than 1 normalized counts are retained.

**To be more specific, where should I apply `filterByExpr` function, before normaliztion (ie on raw counts) or after normaliztion** ?

Code:

group<-factor(paste(targets$Genotype,targets$Time,targets$Treatment,sep="."))
cbind(targets,Group=group)
y<-DGEList(counts = raw_counts, group = group)

#Filterout low count genes
keep<-filterByExpr(y,group = group)
y<-y[keep,keep.lib.size=FALSE]

#Data Normalization
y<-calcNormFactors(y)

#Design matrix
design<-model.matrix(~0+group)
colnames(design)<-levels(group)

#Dsipersion estimation
y<-estimateDisp(y,design)
fit<-glmQLFit(y,design)

Answer: Before normalization (or acutally it does not matter as filterByExp does not consider the normalization since it filters on counts).

What values do you have in the Ctrl1, Ctrl2, Str1 and Str2 coulmns? CPM? In that case the filtering probably works as gene1 have at least 0.657007 * Xe6 reads in two samples. If you have sequenced 40e6 reads that corresponds to at least ~26 reads in two samples.

That said - if you want a more stringent filtering you just have to apply it yourself.
Question: This is a post from [stackoverflow][1] here they show how to extract dedrogram such in form of respective cluster but this is with heatmap.2 function , i m trying to do the same with pheatmap , let's say a group of genes are forming clusters then i want to extract them

Im doing something like this in my code but its giving me all of them togethter not cluster wise as seen in the heatmap

out = pheatmap(data,
color=macolor,

show_rownames = F,cluster_cols=T,cluster_rows=T, scale="row",
cex=1,clustering_distance_rows = "euclidean", cex=1,
clustering_distance_cols = "euclidean", clustering_method = "complete",border_color = FALSE)

res <- data[c(out$tree_row[["order"]]),out$tree_col[["order"]]]

[1]: https://stackoverflow.com/questions/28356354/write-table-with-row-name-ordered-after-clustering-in-rs-heatmap-2

any suggestion or help would be highly appreciated.

Answer: If you add an argument for "kmeans_k" in your pheatmap command, then you can use a very simple option, below:

out = pheatmap(data, color=macolor, kmeans_k = 10, show_rownames = F,cluster_cols=T,cluster_rows=T, scale="row", cex=1,clustering_distance_rows = "euclidean", cex=1, clustering_distance_cols = "euclidean", clustering_method = "complete",border_color = FALSE)

geneclusters <- out[["kmeans"]][["cluster"]]

geneclusters[geneclusters==5]

geneclusters[geneclusters==1]

etc.
Question: Hi there,

I'd be grateful for some help please.

I've got several joint-called trio VCFs (unaffected parents and proband) which I'm analysing. The multi-sample vcf nicely shows 0/0, 0/1, 1/1 for each individual per variant. However, I'd like to convert the notation to three columns of simply het, hom and ref to make it much easier to analyse.

I've checked out the annovar script convert2annovar.pl --withzyg, however this only works for individual sample vcfs which isn't what I need. I'd like all the samples together in the same multi-sample vcf but with zygosity stated for each variant.

Any help or suggestions would be greatly appreciated.
E

Answer: Here is a solution using bioalcidaejdk. http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

$ java -jar src/jvarkit-git/dist/bioalcidaejdk.jar -e 'stream().forEach(V->println(V.getContig()+":"+V.getStart()+" "+V.getGenotypes().stream().map(G->G.getType().name()).collect(Collectors.joining(" "))));' input.vcf

RF01:970 HOM_REF HOM_REF HOM_REF HOM_REF HOM_VAR
RF02:251 HOM_REF HET HET HOM_REF HOM_REF
RF02:578 HOM_REF HOM_REF HOM_REF HOM_VAR HOM_REF
RF02:877 HET HOM_REF HOM_REF HOM_REF HOM_REF
RF02:1726 HOM_REF HET HET HOM_REF HOM_REF
RF02:1962 HET HOM_REF HOM_REF HOM_REF HOM_REF
RF03:1221 HOM_REF HOM_VAR HOM_VAR HOM_REF HOM_REF
RF03:1242 HOM_REF HOM_REF HOM_REF HOM_VAR HOM_REF
RF03:1688 HOM_REF HOM_REF HOM_REF HOM_REF HOM_VAR
RF03:1708 HOM_REF HOM_REF HOM_REF HOM_REF HOM_VAR
Question: I need to align my transcripts with database of plants that I downloaded from UniProt. but it is in txt format and I require FASTA format. PLease suggest a solution. Do I have to retrieve the data from other source?

Thanks

Answer: 1. Go to UniProt website
2. Search for: `taxonomy:viridiplantae`
3. Click on the `Download` button:
- Check `Download all`
- Select **Format:** `FASTA`
- Check `Compressed`
- Click **GO**.
Question: Its a silly question

How to extract list of read names and their corresponding aligned fasta sequence name as table from the reference provided based on alignment.bam file

Answer: Probably I don't understand the question, but what about this?

samtools view alignment.bam | cut -f1,3

Unless fasta identifier != "fasta sequence name"
Question: Hello, everybody!

I am looking for the procedure reversed to cat.

I would like to separate a long fasta-file into many single fasta sequences, so each sequence from the large fa-file will be situated in a new small file.fa each. I usually used something like the following, but it not as good as it looks like... Where is the mistake?

Many thanks for your help!

N.

```
sub get_fasta {
my $inputName = $_[0];
my $fh;
# new input record separator for fasta-files
$/ = "\n>";
open ($fh, $inputName) or die "cannot open file: $!\n";
my @entries = <$fh>;
# remove last character '>' except for the last entry
chop(@entries[0..$#entries - 1]);
# add '>' except for the 1st entry
foreach (@entries[1..$#entries]) {$_ = '>'.$_;}
# print"@entries";
close($fh);
$/="\n";
return \@entries;
```

Answer: linearize your fasta sequence (two lines per sequence) and use `split -l 2`

cat *.fasta |\
awk '/^>/ {if(N>0) printf("\n"); printf("%s\n",$0);++N;next;} { printf("%s",$0);} END {printf("\n");}' |\
split -l 2 - seq_
Question: I am working on a RNA-Seq data. The Fastqc shows that data has duplication level of >65%. Can anyone tell me how I can extract the duplicated reads and quantify them? I am also interested to rank those duplicated reads.

Thanks!

Answer: After marking duplicates (with your favorite program) you can use the samtools flag to pull out the reads you want.

Just set the include flag: samtools view -f 0x400:

0x400 PCR or optical duplicate
Question: I could not find any information on 1) how are built the annotation files that one can find in UCSC database and 2) what's the difference between the different annotation files they offer.

Do someone know if they have their own pipeline to build an annotation ? Can we find instructions of how they're built ?

For example, when trying to download chicken annotation data, I go to https://hgdownload.soe.ucsc.edu/goldenPath/galGal6/bigZips/ and then click "genes/" folder. Here I have the choice between 3 annotations:

- galGal6.ensGene.gtf.gz
- galGal6.ncbiRefSeq.gtf.gz
- galGal6.refGene.gtf.gz

I am not sure which one I should use. And they are very different, even just the features numbers (number of lines in each file) is highly variable:

- 833601 galGal6.ensGene.gtf
- 1768359 galGal6.ncbiRefSeq.gtf
- 163989 galGal6.refGene.gtf

Any advice on that would be welcome.

Answer: Hi NathL,

Those files represent three annotation tables that make-up similarly named UCSC Genome Browser tracks, converted directly into GTF with the `genePredToGtf` command. These files were created since there is an error in Table Browser where the geneId and transcriptId are identical, producing improper GTF files.

Ensembl and NCBI both release independent gene annotations with different quality criteria. The annotations for the refGene (UCSC RefSeq) track were generated from UCSC's realignment of RefSeq RNAs with NM and NR accessions, ignoring the XM and XR unvalidated gene prediction category of annotations. You can look at these specific datasets on the Genome Browser visualization. The choice of gene set depends on your application and tolerance of false positives versus potentially missing annotations. You can read more about their differences in the FAQ pages or the track description pages:

- https://genome.ucsc.edu/FAQ/FAQdownloads.html#download37
- http://genome.ucsc.edu/FAQ/FAQgenes.html#ens
- http://genome.ucsc.edu/cgi-bin/hgTrackUi?db=galGal6&g=refSeqComposite
- http://genome.ucsc.edu/cgi-bin/hgTrackUi?db=galGal6&g=ensGene

Hope that helps!

Daniel Schmelter

UCSC Genome Browser Support Team