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.
•
written
by
Luyi Tian •
You'll want to make sure that the coordinates are 0-based for BED files, so probably subtract 1 from each start and end coordinate.
You have two operations that need to be performed. 1) converting your coordinates into valid input for a tool - either 1-based chr:start-end or 0-based chr start end and 2) transposing the FASTA output so that each entry is on one line. You can do this using pyfaidx:
$ pip install pyfaidx
$ awk '{print $1,$2-1, $2}' records.txt | xargs faidx hg19.fa --bed - --transform transposed
In R/Bioconductor you can use the BSGenome object. For example:
> library("BSgenome.Hsapiens.UCSC.hg19")
> getSeq(BSgenome.Hsapiens.UCSC.hg19, 'chr1', start=10000, end=10200)
201-letter "DNAString" instance
seq: NTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAA...CTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACCCTAACCCTAA
You can combine this with the TxDb object to get the sequences of all human genes:
> library(Homo.sapiens)
> getSeq(BSgenome.Hsapiens.UCSC.hg19, genes(TxDb.Hsapiens.UCSC.hg19.knownGene)[1:10])
A DNAStringSet instance of length 10
width seq names
[1] 16043 ATAGTTACGCAGGCGCAGTGGGGAGAAACGCGACGCCTTGGGCCGCTCTGCCGAATGCAACC...GGCTCTGTTTAAAATGTCTTTGGACTCCCAGGGCTGAGTGGGCTGGGATCTCGTGTCCTCAA 1
[2] 9969 TGAGATCACTTCCCTTGCAGACTTTGGAAGGGAGAGCACTTTATTACAGACCTTGGAAGCAA...TATTGATTTATTGTTGAATTCATAGTAAATTTTTACTGGTAAATGAATAAAGAATATTGTGG 10
[3] 32214 GGCCCGTTAAGAAGAGCGTGGCCGGCCGCGGCCACCGCTGGCCCCAGGGAAAGCCGAGCGGC...CTGTGGCAACTTGTACTGAAAATCTGGTGCTCAATAAAGAAGCCCATGGCTGGTGGCATGCA 100
[4] 226516 GGGGAGCGCCATCCGCTCCACTTCCACCTCCACATCCTCCACCGGCCAAGGTCCCCGCCGCT...AGAATATAAATGATACACCTCTGACCCCAGCGTTCTGAATAAAATGCTAATTTTGGATCTGG 1000
[5] 355352 CTCAAATACACATCACCAAACAAATTTTCTCTATTATTTGGGTAGGCGTGACTGGTTTTCTT...ACTCGAATTCTCACAAGGAAAAGGCCATTAAAGCTCAAGGTGCATTTCAAACTCCAGGCTAC 10000
[6] 15729 TGTGAAATATGAGTTGGCGAGGAAGATCGACCTATTATTGGCCTAGACCAAGGCGCTATGTA...AATTTGTTCATTAAAATTCTCCCAATAAAGCTTTACAGCCTTCTGCAAAGAAGTCTTGCGCA 100008586
[7] 2784 GTGGGAAGTGCGGCTCCTTTCGCGTCCCCCACCCTCTCGGCTCCGCCTGGCAGCAGCTCCGC...ACATTCAAGTTGGGTATAATGACTTTATTTTCCCAGTAATTCCCTAAATAAAACTATTACAA 100009676
[8] 16428 GAAAGAGAACCTGTAAACGCTCTCGGAATTATGGCGGCGGTGGATATCCGAGGTATACTGTA...ACCTACACACTCACATGTTTACTGGTAGATATGTTTAAAAGCAAAATAAAGGTATTTGTATA 10001
[9] 7704 AGAAATCTCCTCAAGCCAGAGCCTGTGCTGTGAGGGGCTTCGGGACCTTGGGGCAGCTCCTG...GTGCACTTTATTTTGTTAGACCTCAATAAAAAAGTAAAAAAAAAAAACAAAAAAAACCAGAA 10002
[10] 57962 GCGGAGGCCACCCAGAGCTCACAGCCTCCTGCCAGCGCGCTCTCTGTTTCTCTGCAGCCCCG...AAATATATATAAAATTTATTGTATGAAAATTAAGCCTCAATAAACGTGATTATAAAAAACAA 10003
Log in to answer this question.
search this site for samtools faidx