Repository navigation
Expand file tree
/
Copy pathhgdp-tgp-fix-psam.R
More file actions
164 lines (142 loc) · 8.23 KB
/
Copy pathhgdp-tgp-fix-psam.R
File metadata and controls
164 lines (142 loc) · 8.23 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
library(tidyverse)
library(genio)
# this processing is very restricted to the specific file provided
# parsing fails without a larger than default guess max
data <- read_tsv( 'gnomad_meta_updated.tsv', guess_max = 5000 )
# existing psam file, only IDs are useful (everything else is missing)
psam <- read_tsv( 'hgdp-tgp-autosomes.psam' )
# remove the initial hash on the first column name
names( psam )[1] <- sub( '^#', '', names( psam )[1] )
# why don't these match?
nrow( data ) # [1] 4150
nrow( psam ) # [1] 4091
# these two IDs are identical, why???
stopifnot( all( data$s == data$project_meta.sample_id ) )
# remove columns that are all NAs (why???)
dim( data ) # [1] 4150 198
data <- data[, colSums( !is.na( data ) ) > 0 ]
dim( data ) # [1] 4150 151
# quite a few more columns were monomorphic (single value)
data <- data[, apply( data, 2, function ( x ) length( unique( x ) ) ) > 1 ]
dim( data ) # [1] 4150 132
# v2 marks only some TGP samples
#table( data$project_meta.v2_release, data$project_meta.title )
## 1000 Genomes Project Human Genome Diversity Project
## FALSE 1328 948
## TRUE 1874 0
#data %>% select( project_meta.subpop_description, hgdp_tgp_meta.Population, population ) %>% unique %>% filter( hgdp_tgp_meta.Population != population )
## project_meta.subpop_description hgdp_tgp_meta.Population population
## 1 Mbuti Pygmy MbutiPygmy Mbuti
## 2 Biaka Pygmy BiakaPygmy Biaka
## 3 Melanesian Melanesian Bougainville
## 4 Papuan Papuan PapuanSepik
## 5 Papuan Papuan PapuanHighlands
## 6 Italian Italian BergamoItalian
## 7 Mongola Mongola Mongolian
## 8 Han Han NorthernHan
#data %>% select( population, longitude, hgdp_tgp_meta.Longitude ) %>% unique %>% filter( longitude != hgdp_tgp_meta.Longitude )
## population longitude hgdp_tgp_meta.Longitude
## 1 PapuanHighlands 145. 143
## 2 Mongolian 119 111
## 3 NorthernHan 108. 114
# manually pruned further, many numbers rounded to the same
#data %>% select( population, latitude, hgdp_tgp_meta.Latitude ) %>% unique %>% filter( latitude != hgdp_tgp_meta.Latitude )
## population latitude hgdp_tgp_meta.Latitude
## 8 PapuanHighlands -6.1 -4
## 14 Mongolian 48.5 45
## 17 NorthernHan 34.7 32.3
# the paper uses `population`, `longitude`, and `latitude`
# the paper uses the last one, and it's the most logical and convenient grouping I'd say
#data %>% select( project_meta.project_pop, hgdp_tgp_meta.Study.region, hgdp_tgp_meta.Genetic.region ) %>% unique
## project_meta.project_pop hgdp_tgp_meta.Study.region hgdp_tgp_meta.Genetic.r…¹
## 1 sas Central_South_Asia CSA
## 2 afr Africa AFR
## 3 oth Oceania OCE
## 4 nfe Europe EUR
## 5 mid Middle_East MID
## 6 amr America AMR
## 7 eas East_Asia EAS
## 8 nfe EUR EUR
## 9 fin EUR EUR
## 10 eas EAS EAS
## 11 amr AMR AMR
## 12 sas SAS CSA
## 13 afr AFR AFR
## # want to preserve their relatedness filters, this fishes out columns describing relatedness
## data %>% select( names( data ) %>% grep( 'rela', x = ., value = TRUE ) ) %>% unique %>% View()
## # the disagreement is in a single pair, columns trade on which one is kept
## data %>% filter( s %in% psam$IID, sample_filters.release_related != sample_filters.all_samples_related ) %>% View
## # count cases; the last column detects a lot more related people, disagrees with the paper counts more
## sum( data$sample_filters.release_related, na.rm = TRUE ) #[1] 711
## sum( !data$sample_filters.release_related, na.rm = TRUE ) #[1] 3408
## sum( data$sample_filters.all_samples_related, na.rm = TRUE ) #[1] 711
## sum( !data$sample_filters.all_samples_related, na.rm = TRUE ) #[1] 3408
## sum( !is.na( data$relatedness_inference.relationships ) ) #[1] 1925
## sum( is.na( data$relatedness_inference.relationships ) ) #[1] 2225
# after a lot of inspection, keep these columns, rename and reorder them to match my old pops-annot.txt files
data <- data %>% select( id = s, sex = project_meta.sex, rel = sample_filters.release_related, pop = population, pop_human = project_meta.subpop_description, superpop = hgdp_tgp_meta.Genetic.region, project = hgdp_tgp_meta.Project, x = longitude, y = latitude )
# TGP includes trios!, HGDP is bigger than my previous WGS set!
table( data$project )
## 1000 Genomes HGDP
## 3202 948
# cd ~/tgp-nygc2
# wc -l *.fam
## 3202 tgp-nygc-autosomes.fam
## 2594 tgp-nygc-autosomes-founders.fam
# cd ../hgdp_wgs/
# wc -l *.fam
## 929 hgdp_wgs_autosomes.fam
# let's compare IDs now that this is a bit visually cleaner
# all individuals with genotypes are in table too
stopifnot( all( psam$IID %in% data$id ) )
indexes <- ! data$id %in% psam$IID
sum( indexes ) # [1] 59
# inspection reveals some random issues with some of these, but unclear in general why they are all excluded, meh
#View( data[ indexes, ] )
# filtered version
data2 <- data[ !indexes, ]
table( data2$project )
## 1000 Genomes HGDP
## 3166 925
table( data2$rel )
## FALSE TRUE
## 3387 704
# paper said it's 3,400 unrelated, but this is closest we get (other column is worse)
# pull the individual data from the pop data
fam <- data2 %>% select( FID = pop, IID = id, SEX = sex, REL = rel )
pop <- data2 %>% select( -id, -sex ) %>% unique
# massage fam data a bit more
# encode sex numerically
fam$SEX <- sex_to_int( substr( fam$SEX, 1, 1 ) )
# sort to match existing psam
fam <- fam[ match( psam$IID, fam$IID ), ]
# add hash to first col name
names( fam )[1] <- paste0( '#', names( fam )[1] )
# first save subset fo unrelated individuals
write_tsv( fam %>% filter( !REL ) %>% select( -REL ), 'unrelated.psam' )
# now save the unfiltered data, just don't include REL (won't need it)
write_tsv( fam %>% select( -REL ), 'hgdp-tgp-autosomes2.psam' )
# now pops
# sort all by superpopulation in this order
superpop_order <- c( 'AFR', 'MID', 'EUR', 'CSA', 'EAS', 'OCE', 'AMR' )
superpop_human <- c(
AFR = 'African',
MID = 'Middle Eastern',
EUR = 'European',
CSA = 'Central South Asian',
EAS = 'East Asian',
OCE = 'Oceanian',
AMR = 'American'
)
pop <- pop[ order( match( pop$superpop, names( superpop_human ) ) ), ]
# add human versions of superpops (only TGP had them before)
pop <- pop %>% mutate( superpop_human = superpop_human[ superpop ] ) %>% relocate( superpop_human, .after = superpop )
# for the nicer human names, HGDP sometimes still has the more obsolete sounding ones
# this calculates the nice names, but puts them in a different table (because we don't want to mess with TGP names)
pop2 <- pop %>% filter( project == 'HGDP', pop != pop_human ) %>% mutate( pop_human = trimws( gsub( "([[:upper:]])", " \\1", pop ) ) )
# this updates the right rows only
indexes <- match( pop$pop, pop2$pop )
indexes2 <- !is.na( indexes )
pop$pop_human[ indexes2 ] <- pop2$pop_human[ indexes[ indexes2 ] ]
# the most painful update still due is the population reordering, which has to be done manually and by trial and error usually, that's not part of this script
write_tsv( pop, 'pops-annot.txt' )