cannin/biostars_qa huggingface dataset:
Question: I have a VCF file with ~30K sites across 131 samples. I am trying to make it include only variant sites, meaning I want to exclude loci where all of my 131 samples have the same genotype, regardless of what the reference allele is. I used GATK SelectVariants with the -env tag, but that only excludes sites where all samples are 0/0, not sites where all samples are 1/1 (homozygous reference.)

I am a pretty terrible coder and struggle to modify VCF files.

My question is: Does anybody have a script or know of a tool that can remove the entire site (line) if all 131 samples (columns?) have 1/1 in the genotype position? Or more generally, if all samples have the same genotype at that site, whether it be 0/0, 0/1, or 1/1 (GATK can do the 0/0 and 0/1, but if it's easier to kill 3 birds with one stone then no problem).

Thanks!

Alex

Answer: using VCFFilterjs: http://lindenb.github.io/jvarkit/VCFFilterJS.html

java -jar dist/vcffilterjs.jar -e 'function accept(v) {var g0= v.getGenotype(0);for(var i=1;i< v.getNSamples();i++) {if(!v.getGenotype(i).sameGenotype(g0)) return true;} return false;}accept(variant);' input.vcf
Question: Hi everyone,

I have two unequal data sets.

A data set
V1 V2 V3
6	42721754	42721769
6	42721757	42721772
6	42721760	42721775
6	42721763	42721778
6	42721766	42721781
6	42721769	42721784
6	42721772	42721787
6	42721775	42721790

B data set
V2	AF
42721757	0.003067485
42721760	0.006134969
42721763	0.006134969
42721766	0.003067485
42721769	0.006134969
42721772	0.006134969
42721775	0.003067485
42721778	0.006134969
42721781	0.003067485
42721784	0.003067485
42721787	0.009202454
42721790	0.009202454

I want to check if the value of `V2` in `B data set` falls between the values of `V2` and `V3` of `A data set`. And when its true, I want value of `AF` in `B data set` to be added as a new column in `A dataset`.
There are two important points.
1. `B data set` is 80 rows while `A data set` is 6000 rows.
2. value of `V2` in `B` may be repeated, which shouldn't get thrown out in final output.

an urgent help with be appreciated.

Thanks,

Answer: Using **R** *data.table*:

library(data.table)

# example data
A <- fread(text = "
V1 V2 V3
6 42721754 42721769
6 42721757 42721772
6 42721760 42721775
6 42721763 42721778
6 42721766 42721781
6 42721769 42721784
6 42721772 42721787
6 42721775 42721790")

B <- fread(text = "
V2 AF
42721757 0.003067485
42721760 0.006134969
42721763 0.006134969
42721766 0.003067485
42721769 0.006134969
42721772 0.006134969
42721775 0.003067485
42721778 0.006134969
42721781 0.003067485
42721784 0.003067485
42721787 0.009202454
42721790 0.009202454")

# set the same keys for both A and B, start and end positions.
# If we have multiple chromosomes, we need to include chrom as key, too.
setDT(A, key = c("V2", "V3"))
setDT(B)
B <- B[, .(V2, V3 = V2, AF)]
setkeyv(B, c("V2", "V3"))

# merge with overlap, check the "type" argument for merge types, see ?foverlaps
res <- foverlaps(B, A)
head(res)
# V1 V2 V3 i.V2 i.V3 AF
# 1: 6 42721754 42721769 42721757 42721757 0.003067485
# 2: 6 42721757 42721772 42721757 42721757 0.003067485
# 3: 6 42721754 42721769 42721760 42721760 0.006134969
# 4: 6 42721757 42721772 42721760 42721760 0.006134969
# 5: 6 42721760 42721775 42721760 42721760 0.006134969
# 6: 6 42721754 42721769 42721763 42721763 0.006134969
Question: Hi Everyone,

How can I extract promoter sequences (ca. 1000bp upstream TSS) in multi-fasta format from genome (also multi-fasta with scaffolds) using information from corresponding GFF file? I've already tried to use GFF-Ex tool, however it didn't help (finished with errors). It is tobacco genome (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/715/135/GCF_000715135.1_Ntab-TN90/).

Does anyone know some other tools for this?

Thanks,

Answer: Finally, I've solved this problem by combining samtools, bedtools as well as custom R script. The pipline placed into bash script is available [here](https://github.com/RimGubaev/extract_promoters).
Question: Hi all,
I am struggling to make a contrast with EdgeR.

design <- model.matrix(~Gender + Age + Group, data = patient.table)

The (disease) Group is a factor with [Control, Mild, Severe]. I want to compare C-M. M-S and C-S. Thus s i

DGEL <- edgeR::DGEList(expressionData)
keep <- edgeR::filterByExpr(DGEL)
DGEL <- DGEL[keep, , keep.lib.sizes=FALSE]
DGEL <- edgeR::calcNormFactors(DGEL, method = "TMM")

DGEL <- edgeR::estimateDisp(DGEL, design)
fit <- edgeR::glmQLFit(DGEL, design)

contrasts <- makeContrasts(c.vs.m = GroupC - GroupM,
m.vs.s = GroupM - GroupS,
c.vs.s = GroupC - GroupS,
levels = design)

However, there is no GroupC in the model matrix. There is only an Intercept, GroupM and GroupS. How can I compare anything against the control group on a full-rank design matrix.

I already searched in the user manual and here on the website, but could not find it. Also, I tried

design <- model.matrix(~Gender + Age + Group,
data = patient.table,
contrasts.arg = list(
Group = contrasts(patient.table$Group, contrasts = FALSE)
)
)

But edgeR did not like that and informed me that the table was not of full rank.
I am pretty sure this is something basic, but I just don't get it as a wet-lab worker.

Answer: Check [edgeR Users Guide][1], chapter 3.3, on how to define the contrasts for such a complex design. For me, the approach that works better (as in, it is easier to understand) is the one described in section 3.3.1. Start by creating a composite `group` variable:

group <- factor( paste0( patient.table$Gender, paste0( patient.table$Age, patient.table$Group ) )

Then specify the design without an intercept, and create the contrasts of interest.

design <- model.matrix( ~0 + group )
colnames( design ) <- levels( group )
contrasts <- makeContrasts( GenderM_Age0_GroupC_vs_M = GenderMAge0GroupC - GenderMAge0GroupM,
# and so on for the other contrasts
levels=design)

To perform the tests:

fit <- glmQLFit( DGEL, design )
qlf <- glmQLFTest( fit, contrast = contrastss[ , "GenderM_Age0_GroupC_vs_M" ] )

This tests for **GroupC** vs **GroupM** differences within **GenderM** and **Age0**.

See https://www.biostars.org/p/453485/ for a recent similar post, you may find the answer there useful.

If you want more complex tests, such as any genes changing over all ages, read sections 3.3.3 and 3.3.4 from the guide.

[1]: http://bioconductor.org/packages/3.12/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf
Question: Hi Guys,
Im new to RNA-seq and R-programming so forgive my ignorance in advance! I'm currently using a programme/script to help me map tRNAs (the supplied notes don't explain in much detail), and after tRNA counts are generated in all my conditions, they use SizeFactors to normalize the dataset (in DESeq2). I've tried to read up on what exactly SizeFactors are and I don't understand it. Could anyone give me an easy to understand definition of what size factors are and why they're used to normalize the data.

Answer: A size factor relates to how many reads there are in each library. One can imagine that if you had two sample where 10% of the reads in each sample were from gene A, but in one sample 1M reads have been sequenced and in the other sample 2M reads had been sequenced then there would be a two fold increase in the counts from gene A in sample 2 compared to sample 1, but the actaul expression levels were the same.

Early RNAseq analysis divided counts by the total number of reads in each library, but this is poor practice for two reasons.

1. Using division means that you lose the discrete nature of the gene counts, and thus negative bionomial statisitcs no longer apply. Thus normalising factors are used as offsets in the linear model, rather than divisors.
2. In most RNAseq samples the most higly expressed genes take up the majority of the reads. Thus in a 1M read sample, if 300k reads came from gene A (leaving 700k for all the others), and that gene doubled in expression to 600k reads (leaving 400k for all the others), the expression of the other genes would appear to half, even though they have stayed the same.

Thus sizeFactors are related to the library size (total number of reads in the library), but are calculated in such a way as to compensate for effect 2 above. One common method (and the one I believe that is used by DESeq2), is to find the 75th centile of the distribution of read counts for each sample, and then calculate a normalisation factor such that the 75th centile is the same across all samples.
Question: Hello everyone,

I want to use awk to count number of values within a range in a specific column, and to print ranges with the count in output. For example, 4th column of my input is like this

1000
4000
5001
7050
4000
3500

With bin size of 5000, I want to know the number of values in range 1 - 5000, 5001 - 10000 and so on. My values are greater than 100,000. I want to count the values as many times as they appear. For example 4000 should be counted twice. The output should be showing two columns showing ranges and the count like

1000-5000 4
5000-10000 2

I found many relevant answers on different forums but nothing exactly serves my purpose. I found one code but it is counting duplicate values only once.

Thanks

Answer: Using **awk**:

$ awk '{print (int($1/5000)+1)*5000}' temp.txt | sort -g | uniq -c
# 4 5000
# 2 10000

Using **R**:

# example input
x <- c(1000, 4000, 5001, 7050, 4000, 3500)

# using cut into 5K intervals
table(cut(x, seq(0, 10000, 5000)))
# (0,5e+03] (5e+03,1e+04]
# 4 2

# using round
table(round(x, -4))
# 0 10000
# 4 2
Question: In this graph 81 cells have being ordered by pseudo time. For example cell s0.76 has been placed at the start of pseudo time and cell s0.11 placed at the end of pseudo time.

![enter image description here][1]

This is the expression profile of this gene on this pseudo time

> Data

s0.76 s0.52 s0.72 s0.36 s0.18 s0.60 s0.54 s0.65 s0.56 s0.4 s0.14 s0.67 s0.1 s0.22 s0.28 s0.80 s0.16 s0.32
6354.316 2376.278 8037.063 4879.946 4989.872 5382.370 13239.227 12085.752 5510.481 9311.501 5303.664 13480.634 9085.402 9585.708 12678.642 8595.338 13901.138 9809.339
s0.8 s0.55 s0.71 s0.41 s0.40 s0.48 s0.68 s0.81 s0.29 s0.78 s0.34 s0.42 s0.13 s0.9 s0.7 s0.24 s0.62 s0.35
12554.518 11001.622 6676.949 9266.956 5303.998 13405.392 8327.440 6903.236 9319.111 4427.627 11029.422 5370.687 4854.454 13199.522 5125.644 11687.103 7546.622 8656.206
s0.23 s0.21 s0.37 s0.25 s0.79 s0.64 s0.44 s0.15 s0.43 s0.63 s0.10 s0.6 s0.19 s0.31 s0.38 s0.39 s0.73 s0.59
4578.953 15771.795 3217.201 6362.859 11523.728 5977.613 4805.703 5373.737 2393.996 9104.929 3532.215 7974.742 4949.271 3285.566 10896.503 3746.348 4296.337 6510.001
s0.2 s0.12 s0.61 s0.66 s0.51 s0.58 s0.50 s0.27 s0.77 s0.69 s0.33 s0.17 s0.74 s0.26 s0.57 s0.49 s0.70 s0.53
3728.896 8761.615 3896.014 2618.151 4106.299 11898.634 2978.410 3129.880 5571.995 4314.655 3867.158 2209.066 1723.784 2741.320 4570.466 5506.473 3120.334 7284.682
s0.30 s0.45 s0.5 s0.47 s0.46 s0.75 s0.20 s0.3 s0.11
9454.997 7446.095 4765.971 7861.282 10393.357 4809.735 12421.501 12615.785 22238.399

Now, my question is, how I could color the plot by knowing a set of cells are in G2 and another set are in S/M phase of cell cycle like this picture?

![enter image description here][2]

These cells are at G2

[1] "s0.1" "s0.2" "s0.3" "s0.6" "s0.7" "s0.8" "s0.9" "s0.10" "s0.12" "s0.13" "s0.14" "s0.15" "s0.16" "s0.19" "s0.21" "s0.22" "s0.23" "s0.24" "s0.25" "s0.28" "s0.29" "s0.31"
[23] "s0.32" "s0.34" "s0.35" "s0.37" "s0.38" "s0.39" "s0.42" "s0.43" "s0.44" "s0.48" "s0.50" "s0.51" "s0.54" "s0.55" "s0.56" "s0.59" "s0.60" "s0.61" "s0.62" "s0.63" "s0.64" "s0.65"
[45] "s0.66" "s0.67" "s0.68" "s0.70" "s0.71" "s0.73" "s0.78" "s0.79" "s0.80" "s0.81"

these cells are in M/S

[1] "s0.4" "s0.5" "s0.11" "s0.17" "s0.18" "s0.20" "s0.26" "s0.27" "s0.30" "s0.33" "s0.36" "s0.40" "s0.41" "s0.45" "s0.46" "s0.47" "s0.49" "s0.52" "s0.53" "s0.57" "s0.58" "s0.69"
[23] "s0.72" "s0.74" "s0.75" "s0.76" "s0.77"

I tried to do that by this codes but the order of cells is not based on the pseudo time

> head(data)
value cell capture
1 30.59486 s0.76 G2
2 78.12876 s0.52 G2
3 63.36236 s0.72 G2
4 54.48303 s0.36 m/s
5 196.49463 s0.18 m/s
6 59.76789 s0.60 G2
>

![enter image description here][3]

[1]: https://image.ibb.co/gEKXsK/Rplot.jpg
[2]: https://image.ibb.co/c5PJzz/btw372f3p.jpg
[3]: https://image.ibb.co/gptvCK/Rplot01.jpg

> p=ggplot(data, aes(x=cell, y=value, color=capture))
> p+geom_point(alpha=.3)

Answer: p+geom_point(alpha=.3)+scale_x_discrete(limits=as.factor(data$cell))
Question: Recently, I searching paper about time series or time course, I found that these two words equally used by researchers. However, time series seems used more wider. wikipedia also use time series (https://en.wikipedia.org/wiki/Time_series). What the difference between these two phrases? Is there any specific context that one phrases more frequently used?

Best,

Pengcheng

Answer: They're the same thing. You'll find that there are many words and phrases in English that mean the same thing but happen to be used more in one group than another (e.g., in my original field of study, a "time series analysis" could sometimes be called a "non-stationary analysis").
Question: We are in the process of developing analysis pipelines for WGS,RNA and Methyl-seq data. The first projects considers 100 patients with a common disease. All the analysis would be run on data from these patients. Our institute's computational center is offering us 5 nodes (24 CPU core, 2.6GHz clock speed, 128GB RAM memory each), I believe that necessary storage space would also be provided to us. Based on your experiences would you please tell me if having 5 nodes is going to be sufficient, I understand that there cannot be one single answer to this question but any help from you all is much appreciated.

Answer: The two big things are

1) How much ram does a single instance of the "largest" program you'll run need?

This is the only hard requirement. When I've run trinity I've needed 1TB of ram before. If I didn't have a machine with 1TB of RAM, I couldn't have run Trinity. If the program with the highest memory consumption will fit on 124gb, then that's fine.

2) How long are you willing to wait for things to run?

More nodes means you can run more stuff at once. You don't really *need* them. Figure out which is more important to you, time or money? Spend more money, get more nodes, higher throughput. Save money, get fewer nodes, spend more time waiting for results.
Question: I have f3 statistics in the form f3(Outgroup:POI_1, Comparative) and f3(Outgroup:POI_2, Comparative) and so on.

For a handful of pairs of POIs (populations of interest) and over a hundred comparative populations

For each pair of POIs, I want to make a scatter plot, where one f3 value is on the X axis and the others is on the Y axis, but I want both sets of standard errors graphed.

Basically, I want some conceptually similar to this plot from Allentoft et al. 2015 (http://www.nature.com/nature/journal/v522/n7555/fig_tab/nature14507_SF2.html).

Answer: This is really an R "question", but here you go.
Question: Hello biostars,

I am new to this type of analysis. I have IDR.NarrowPeak data and I would like to normalise it to FPKM using R.
I read that the packages limma and DESeq are helpful but I still have problems creating the count matrix.
Is there a step I am missing?

Any help/suggestion will be appreciated.

Answer: You should try DiffBind, which uses DESeq and edgeR in the background and automates this process: https://bioconductor.org/packages/release/bioc/html/DiffBind.html
Question: hi there,
is there a way to get count of SNP, indels, CNVs etc from a VCF file, so some thing like

SNPs = ?

Insertions = ?

Deletions = ?

CNVs = ?

using simple linux commands

thanks,
a

Answer: There are a couple of ways that variant type is annotated within a VCF file, so there are correspondingly a few ways to get close to what you want. Here's one choice that should work with most VCF files:

Use the `vcftools` tool `vcf-annotate` to fill in the variant type field:

zcat in.vcf.gz | vcftools_0.1.9/bin/vcf-annotate --fill-type > out.vcf

Then count up the variants by looking at the (newly-filled) TYPE field:

grep -oP "TYPE=\w+" out.vcf | sort | uniq -c

Or in one step that doesn't change the original VCF file:

zcat in.vcf.gz | vcftools_0.1.9/bin/vcf-annotate --fill-type | grep -oP "TYPE=\w+" | sort | uniq -c

On an example I had, this yielded:

3410 TYPE=del
4487 TYPE=ins
56744 TYPE=snp

1000 Genomes VCF files will be annotated in a finer-grained way (e.g. choices including DUP, INV, CNV, TANDEM, see [here](http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-41)), but I'm not sure how to get their range of annotations from your own raw read data. However, if these distinctions are critical to you, that may be a useful direction to explore.
Question: Dear All,

I have found a total list of about 3000 transcripts that I am interested in their gene ontology functions. I have the .fasta file of each transcripts. Can anyone advise on how should I go about to extract the GO symbols for each genes?

The pipeline I can think of is input the .fasta file onto GUI-version of Blast2Go and run a remote blast, and then map the GO from Blast2Go. I am already running this right now but it seems quite slow.

Any alternative strategies here?

regards

Ziyi

Answer: Blast2GO is impractially slow in my experience and doesn't fit well into a Unix pipeline, though that is not the usual target usage. Another approach is HMMER2GO if you don't find a few commands. Here is an example starting with some transcript assemblies called "genes.fasta" and ending with GO terms for each gene/contig.

hmmer2go getorf -i genes.fasta -o genes_orfs.faa
hmmer2go run -i genes_orfs.faa -d Pfam-A.hmm
hmmer2go mapterms -i genes_orfs_Pfam-A.tblout -o genes_orfs_Pfam-A_GO.tsv --map

There is some documenation on the wiki or at the command line and a full example, as well.
Question: I'm working from the Biostar Handbook and trying to do the alignment with tophat. I'm on page 480 for reference.

I was able to successfully run

$ bowtie2-build $REF $IDX

However, then it says that I need to change the invocation of the aligner to

$ tophat -G $GTF -o tophat_hbr1 $IDX $R1 $R2
bash: /usr/bin/tophat: No such file or directory
(bioinfo)

And that's my output. If I do which tophat, nothing comes up, so it seems like it's just not there. So, I tried to install tophat with

$ conda install tophat
Fetching package metadata .................
Solving package specifications: .

UnsatisfiableError: The following specifications were found to be in conflict:
- python 3.6*
- tophat -> python 2.7*
Use "conda info " to see the dependencies for each package.

(bioinfo)

I saw some solutions online for changing a line of code in tophat to do this, but I don't even know where that file is. I can downgrade python, but I'm not sure if that will work, and I'm worried that I'll break other things if I do. Also, when I do

$ ls /usr/bin | grep python
dh_python2
dh_python3
python
python-config
python2
python2-config
python2.7
python2.7-config
python3
python3.4
python3.4m
python3m
x86_64-linux-gnu-python-config
x86_64-linux-gnu-python2.7-config
(bioinfo)
moltres@moltres-ao ~/biostar/Sequencing/griffith
$ python -V
Python 3.6.3
(bioinfo)

I see a bunch of different python versions, and I'm not sure if I should switch between them or what.

I'm a little frustrated because I bought this book because I'm a beginner, but I keep running into errors that are over my head. I'd really appreciate help here.

Answer: [Create a new environment for Tophat][1]:

conda create -n tophat2 tophat2

should create a new environment and install Tophat2 and its correct dependencies.

[1]: https://conda.io/docs/user-guide/tasks/manage-environments.html#creating-an-environment-with-commands
Question: Hello,

I'm new here and and I need help please!

I need to download a list of all human genes with their respective gene symbol | chromosome | strand | transcription start site | Txen | and **Ensembl gene name**,

Actually I'm using UCSC table and I get something like this:

```
#hg38.knownGene.name hg38.knownGene.chrom hg38.knownGene.strand hg38.knownGene.txStart hg38.knownGene.txEnd hg38.kgXref.geneSymbol
uc001aaa.3 chr1 + 11873 14409 DDX11L1
```

And I want to know if it's possible to include ensembl gene symbol with UCSC table or with another method

Thank you in advance

Cherif

Answer: Unfortunately UCSC table browser doesn't have Ensembl gene track for hg38. But they do have it for hg19 and the below command should work for you.

```
mysql \
--user=genome \
-N \
--host=genome-mysql.cse.ucsc.edu \
-A \
-D hg19 \
-e "select ensGene.name, name2, chrom, strand, txStart, txEnd, value from ensGene, ensemblToGeneName where ensGene.name = ensemblToGeneName.name" > \
output.txt
```

You can try the same command with hg38 but you will have to choose other gene models such as refseq or ucsc.
Question: Hi all, I have several millions of records with format `chromosome position` like:

```
chr1 10001
chr1 144244
chr2 ...
```

And I know it is from human hg19. Now I want to get the ATGC info for each position so it would become:

```
chr1 10001 A
chr1 144244 G
chr2 ...
```

Anyone know what is the best practice for such job. Tools like UCSC table work for single record, but too slow for a lot of them.

Answer: `bedtools getfasta` will also work for this job. You should create a bed file with the chromosome start and end positions that you'd like to extract:

chr1 10001 10002
chr1 144244 144245
...
Question: There are so many things to read about ga4gh: I'm a little bit lost. So here are a few questions about [GA4GH][1] / Beacons :

- I've got some VCF files. What is the correct way to expose my data? How can I secure my data for my collaborators? Is there a java server implementation?
- If I build a beacon server from scratch (learning... learning...) is there a tool that can be used to test if my implementation is [ga4gh][1] compliant?
- Where should I declare that I'm providing a beacon server?
- The schema use [avro][2]; Is it possible to return plain JSON or does using Avro means that I need to serialize the data as a binary stream? Is there any API converting jaav+htsjdk to avro?
- What's the best place to ask questions about ga4gh?

[1]: http://genomicsandhealth.org/
[2]: https://avro.apache.org/

Answer: This should be the easiest beacon server to run (disclosure: I'm the developer): https://github.com/maximilianh/ucscBeacon

Should run on any OS, brings its own webserver or runs as a CGI. Can import directly: VCF, Complete Genomics, BED, HGMD and LOVD formats. A single python script, no dependencies. Let me know if something doesn't work for you.
Question: Hi,

I have to plot allele frequencies of two different SNP chip datasets. I have two VCF files and would like to make a scatterplot in which these 2 datasets are plotted one against each other. What is the easiest way to do this?
I apologize if this question was already posted.

Answer: bcftools query -i 'N_ALT=1' -f '%CHROM\_%POS\_%REF\_%ALT\t%INFO/AF\n' in1.vcf.gz | sort -t $'\t' -k1,1 > a.tsv
bcftools query -i 'N_ALT=1' -f '%CHROM\_%POS\_%REF\_%ALT\t%INFO/AF\n' in2.vcf.gz | sort -t $'\t' -k1,1 > b.tsv
join -a 1 -a 2 -o '1.2,2.2' -t $'\t' -1 1 -2 1 -e 0.0 a.tsv b.tsv > compare.tsv

plot with R : https://stackoverflow.com/questions/42232340/how-to-read-tab-separated-files-in-r-to-generate-scatterplot
Question: I'm a Computer Science undergraduate writing my thesis which is about analyzing open source bioinformatics projects and extracting standalone modules. I need those projects to be written in Java and preferably be stable or mature. So far I've found a bunch of projects via SourceForge and I was wondering if there are any notable ones that are not hosted there. A different repository that has a bioinformatics category would be helpful too.

Here's a list of the most promising projects I have so far:

- Jenetics
- JGap
- Jmol
- Juicebox
- Picard
- RDP Classifier
- Shim
- Tassel
- The Chemistry Development Kit
- Toxfree
- TreeView
- VarScan

Thanks in advance.

Answer: - [Archaeopteryx][1]
- [Cytoscape][2]

[1]: https://sites.google.com/site/cmzmasek/home/software/archaeopteryx
[2]: http://www.cytoscape.org/