64ef1345bfa1476ad1c2a1ad7cf413e95556966b angie Wed Mar 17 11:53:07 2021 -0700 We can actually end up with more bioSamples than we're looking for due to including past results, so compare IDs instead of counts. diff --git src/hg/utils/otto/sarscov2phylo/getNcbi.sh src/hg/utils/otto/sarscov2phylo/getNcbi.sh index cf77efe..96e63f0 100755 --- src/hg/utils/otto/sarscov2phylo/getNcbi.sh +++ src/hg/utils/otto/sarscov2phylo/getNcbi.sh @@ -26,55 +26,53 @@ > ncbi_dataset.zip rm -rf ncbi_dataset unzip ncbi_dataset.zip # Creates ./ncbi_dataset/ # This makes something just like ncbi.datasets.tsv from the /table/ API query: jq -c -r '[.accession, .biosample, .isolate.collectionDate, .location.geographicLocation, .host.sciName, .isolate.name, .completeness, (.length|tostring)] | join("\t")' \ ncbi_dataset/data/data_report.jsonl \ | sed -e 's/COMPLETE/complete/; s/PARTIAL/partial/;' \ | sort \ > ncbi_dataset.tsv # Use EUtils (esearch) to get all SARS-CoV-2 BioSample GI# IDs: $scriptDir/searchAllSarsCov2BioSample.sh +sort all.biosample.gids.txt > all.biosample.gids.sorted.txt # Copy yesterday's all.bioSample.txt so we don't have to refetch all the old stuff. if [ -e ../ncbi.$prevDate/all.bioSample.txt.xz ]; then xzcat ../ncbi.$prevDate/all.bioSample.txt.xz > all.bioSample.txt - grep ^Accession all.bioSample.txt | sed -re 's/^.*ID: //' | sort > ids.loaded - sort all.biosample.gids.txt > all.biosample.gids.sorted.txt + grep ^Accession all.bioSample.txt | sed -re 's/^.*ID: //' | sort -u > ids.loaded comm -23 all.biosample.gids.sorted.txt ids.loaded > ids.notLoaded else cp -p all.biosample.gids.txt ids.notLoaded fi wc -l ids.notLoaded # Use EUtils (efetch) to get BioSample records for the GI# IDs that we don't have yet: time $scriptDir/bioSampleIdToText.sh < ids.notLoaded >> all.bioSample.txt -gidCount=$(wc -l < all.biosample.gids.txt) -accCount=$(grep ^Accession all.bioSample.txt | sort -u | wc -l) -if [ $gidCount != $accCount ]; then - echo "Number of Accession lines ($accCount) does not match number of numeric IDs ($gidCount)" + grep ^Accession all.bioSample.txt | sed -re 's/^.*ID: //' | sort > ids.loaded - sort all.biosample.gids.txt > all.biosample.gids.sorted.txt comm -23 all.biosample.gids.sorted.txt ids.loaded > ids.notLoaded +if [ -s ids.notLoaded ]; then echo Retrying queries for `wc -l < ids.notLoaded` IDs $scriptDir/bioSampleIdToText.sh < ids.notLoaded >> all.bioSample.txt - accCount=$(grep ^Accession all.bioSample.txt | sort -u | wc -l) - if [ $gidCount != $accCount ]; then + grep ^Accession all.bioSample.txt | sed -re 's/^.*ID: //' | sort > ids.loaded + comm -23 all.biosample.gids.sorted.txt ids.loaded > ids.notLoaded + if [ -s ids.notLoaded ]; then echo "Still have only $accCount accession lines after retrying; quitting." exit 1 fi fi # Extract properties of interest into tab-sep text: $scriptDir/bioSampleTextToTab.pl < all.bioSample.txt > all.bioSample.tab # Extract BioSample tab-sep lines just for BioSample accessions included in the ncbi_dataset data: tawk '$2 != "" {print $2;}' ncbi_dataset.tsv \ | grep -Fwf - all.bioSample.tab \ > gb.bioSample.tab # Use BioSample metadata to fill in missing pieces of GenBank metadata and report conflicting # sample collection dates: