Get the gene2pubmed database day by day

python scripts/pmids_by_date.py --startdate 1990/01/01 --enddate 2017/11/05; \

Consolidate the publications from each day into one complete list

rm ~/data/genbank-data/hg19/recent_pmid_year.ssv; for file in $(find data/pmid_by_date -name "*.ssv"); do cat $file | awk '{split($1,a,"-"); print a[1], $2}' >> ~/data/genbank-data/hg19/recent_pmid_year.ssv; done;

get the gene citation table

wget -N -P ~/data/genbank-data/ ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2pubmed.gz; gunzip ~/data/genbank-data/gene2pubmed.gz; mkdir ~/data/genbank-data/hg19/; \

filter out human genes

zcat ~/data/genbank-data/gene2pubmed.gz | awk '{if ($1 == 9606) print;}' > ~/data/genbank-data/hg19/gene2pubmed cat ~/data/genbank-data/gene2pubmed | python scripts/genes_by_species.py - | sort -nk 2 > results/genes_by_species.tsv; cp results/genes_by_species.tsv ~/Dropbox/tmp/graphic-idea-3.tsv

get chromosome lists

wget -N -P ~/data/ucsc-data/hg19/ http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/chromInfo.txt.gz

get gene exon and intron information

wget -N -P ~/data/ucsc-data/hg19/ http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/refGene.txt.gz

get the gene information table

wget -N -P ~/data/genbank-data/ ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene_info.gz; gunzip ~/data/genbank-data/gene_info.gz; mkdir ~/data/genbank-data/hg19/; \

extract human data

zcat ~/data/genbank-data/gene_info.gz | awk '{if ($1 == 9606) print;}' > ~/data/genbank-data/hg19/gene_info

get gene to refseq information so that we can use it to get transcript level information

wget -N -P ~/data/genbank-data/ ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2refseq.gz; mkdir ~/data/genbank-data/hg19/; zcat gene2refseq.gz | awk '{if ($1 == 9606) print;}' > ~/data/genbank-data/hg19/gene2refseq

# aggregate by citation counts and by year
python scripts/summarize_gene_citations.py; \
python scripts/summarize_gene_citations_all_species.py; \
cp gene_info_by_year.tsv ~/Dropbox/tmp/graphic-idea-2.tsv; \

# get a ranking of human genes
# created using a different pipeline than the one in summarize_gene_citations
# and described here
# https://hms-dbmi.github.io/higlass-docs/gene_annotations.html#creating-the-track
# this is a good file to cross-reference with the the ones above to check that the
# values are correct
cp ~/data/genbank-data/hg19/geneAnnotationsExonUnions.bed ~/Dropbox/tmp/graphic-idea-1.tsv \

    awk '{if ($4 == "TP53") print}' gene_info_by_year.tsv | sort -nk 7

calculate how many of the papers referencing human genes reference the top 100 genes

get the top 100 genes

sort -nk 5 ~/data/genbank-data/hg19/geneAnnotationsExonUnions.bed | tail -n 100 | awk '{print $8}' > ~/data/genbank-data/hg19/top100_genes

and for the top 10 genes

sort -nk 5 ~/data/genbank-data/hg19/geneAnnotationsExonUnions.bed | tail -n 10 | awk '{print $8}' > ~/data/genbank-data/hg19/top10_genes

#532451 python scripts/pubs_for_genes.py ~/data/genbank-data/gene2pubmed ~/data/genbank-data/hg19/top100_genes python scripts/pubs_for_genes.py ~/data/genbank-data/gene2pubmed ~/data/genbank-data/hg19/top10_genes

####### cat ~/data/genbank-data/${ASSEMBLY}/gene2pubmed | awk '{print $2}' | sort | uniq -c | awk '{print $2 "\t" $1}' | sort > ~/data/genbank-data/${ASSEMBLY}/gene2pubmed-count head ~/data/genbank-data/${ASSEMBLY}/gene2pubmed-count

output -> geneid \t refseq_id

cat ~/data/genbank-data/${ASSEMBLY}/gene2refseq | awk -F $'\t' '{ split($4,a,"."); if (a[1] != "-") print $2 "\t" a[1];}' | sort | uniq > ~/data/genbank-data/${ASSEMBLY}/geneid_refseqid head ~/data/genbank-data/${ASSEMBLY}/geneid_refseqid wc -l ~/data/genbank-data/${ASSEMBLY}/geneid_refseqid

#output -> geneid \t refseq_id \t citation_count

join ~/data/genbank-data/${ASSEMBLY}/geneid_refseqid ~/data/genbank-data/${ASSEMBLY}/gene2pubmed-count | sort -k2 > ~/data/genbank-data/${ASSEMBLY}/geneid_refseqid_count

head ~/data/genbank-data/${ASSEMBLY}/geneid_refseqid_count wc -l ~/data/genbank-data/${ASSEMBLY}/geneid_refseqid_count

output -> geneid \t refseq_id \t chr (5) \t strand(6) \t txStart(7) \t txEnd(8) \t cdsStart(9) \t cdsEnd (10) \t exonCount(11) \t exonStarts(12) \t exonEnds(13)

join -1 2 -2 2 ~/data/genbank-data/${ASSEMBLY}/geneid_refseqid_count ~/data/genbank-data/${ASSEMBLY}/sorted_refGene | awk '{ print $2 "\t" $1 "\t" $5 "\t" $6 "\t" $7 "\t" $8 "\t" $9 "\t" $10 "\t" $11 "\t" $12 "\t" $13 "\t" $3; }' | sort -k1 > ~/data/genbank-data/${ASSEMBLY}/geneid_refGene_count

head ~/data/genbank-data/${ASSEMBLY}/geneid_refGene_count wc -l ~/data/genbank-data/${ASSEMBLY}/geneid_refGene_count

output -> geneid \t symbol \t gene_type \t name \t citation_count

join -1 2 -2 1 -t $'\t' ~/data/genbank-data/${ASSEMBLY}/gene_info ~/data/genbank-data/${ASSEMBLY}/gene2pubmed-count | awk -F $'\t' '{print $1 "\t" $3 "\t" $10 "\t" $12 "\t" $16}' | sort -k1 > ~/data/genbank-data/${ASSEMBLY}/gene_subinfo_citation_count head ~/data/genbank-data/${ASSEMBLY}/gene_subinfo_citation_count wc -l ~/data/genbank-data/${ASSEMBLY}/gene_subinfo_citation_count

1: chr (chr1)

2: txStart (52301201) [9]

3: txEnd (52317145) [10]

4: geneName (ACVRL1) [2]

5: citationCount (123) [16]

6: strand (+) [8]

7: refseqId (NM_000020)

8: geneId (94) [1]

9: geneType (protein-coding)

10: geneDesc (activin A receptor type II-like 1)

11: cdsStart (52306258)

12: cdsEnd (52314677)

14: exonStarts (52301201,52306253,52306882,52307342,52307757,52308222,52309008,52309819,52312768,52314542,)

15: exonEnds (52301479,52306319,52307134,52307554,52307857,52308369,52309284,52310017,52312899,52317145,)

join -t $'\t' ~/data/genbank-data/${ASSEMBLY}/gene_subinfo_citation_count ~/data/genbank-data/${ASSEMBLY}/geneid_refGene_count | awk -F $'\t' '{print $7 "\t" $9 "\t" $10 "\t" $2 "\t" $16 "\t" $8 "\t" $6 "\t" $1 "\t" $3 "\t" $4 "\t" $11 "\t" $12 "\t" $14 "\t" $15}' > ~/data/genbank-data/${ASSEMBLY}/geneAnnotations.bed head ~/data/genbank-data/${ASSEMBLY}/geneAnnotations.bed wc -l ~/data/genbank-data/${ASSEMBLY}/geneAnnotations.bed

python scripts/exonU.py ~/data/genbank-data/${ASSEMBLY}/geneAnnotations.bed > ~/data/genbank-data/${ASSEMBLY}/geneAnnotationsExonUnions.bed wc -l ~/data/genbank-data/${ASSEMBLY}/geneAnnotationsExonUnions.bed