편집 및 솔루션
내 원래 질문의 표현이 형편없었고 바퀴를 재발명하려고 했기 때문에 이제 내 질문에 답하고 있습니다(다른 사람에게 도움이 될 수도 있음).
gff2fasta는 제가 필요한 작업을 정확하게 수행하는 도구입니다. 즉, 완전한 게놈(FULLGENOME.fasta라는 fasta 형식의 거대한 파일)에서 주어진 DNA 단편을 추출하는 것입니다.
원하는 조각이 어디에 있는지 안다면 TEST.gff라는 파일을 생성하여 gff2fasta 조각(여기: sca_5_chr5_3_0)과 시작(여기: 2390621) 및 끝(여기: 2391041)의 스캐폴딩을 지정합니다. 예를 들어:
sca_5_chr5_3_0 JGI CDS 2390621 2391041 . + 0 name "e_gw1.5.88.1"; proteinId 40463;
그럼 그냥 달려야지
gff2fasta.pl -gff TEST.gff -fasta FULLGENOME.fasta -out test
그런 다음 gff2fasta는 TEST.gff에서 정보를 가져와 다음과 같은 test라는 파일에 원하는 비트를 출력합니다.
>sca_5_chr5_3_0_2390621_2391041_F
ATGACGACCCGCGCACCCAAAGACACATACGCTCAGCCCGACTATGAGGAGGCTCACCTAGCGACGTTTGCAGCCCCAAA
AGGCTACCCTATCGAGTCTATGCTACCCCCTAGCGTGAAGAGGGAGACCTTTGAACAGGCCCTAGCCGAGTTTACCGACG
CTATCGGCAAAGATTATGTCTTTATCGGCGATGCTCTCTCTCACTACATCGATCCCTACGACATCTATATCGATGATAGT
GAGAGAAGGAGGATGCCGAGCGCGGCTGTTTGGTACGTAACGCGGCACCCATCGAAACCTAGCAGCACTGACAAGTTTCC
GCGTCTAGTCCCTCTTCGCTCGAGGAAGCTAAGCAGGCTCTCAAGGTTGCGAACAAATACGGCATTCCGATCTGGGCATT
TTCCAGAGGCAAGAACCTGGG
도움을 주신 terdon과 다른 분들께 감사드립니다!
-------------
이것은 더 많은 정보와 세부 사항을 포함하는 해당 질문의 연속입니다.unix: 파일에서 문자 10~80 가져오기
거의 다 왔다고 생각하지만 아직 도움이 필요합니다.
최대한 명확하게 설명하려고 노력하고 있지만, 이는 매우 기술적인 질문이라는 것을 알고 있으므로 더 명확하게 설명할 수 있는 부분이 있으면 알려주세요!
내가 가지고 있는 파일은 세 개입니다.
관리자 참고사항: 파일을 업로드할 수 있나요? 어떻게 해야할지 모르겠어...
하나문서(N_haematococca_1.fasta) 다음에서 이름을 추출해야 합니다.
head -1 N_haematococca_1.fasta | cut -f4 -d'|'
이 경우 이름은 다음과 같습니다.
e_gw1.5.88.1
질문 1: 위의 코드는 잘 작동하지만 나중에 사용할 수 있는 변수에 이름(e_gw1.5.88.1)을 저장하는 데 문제가 있습니다...
해당 이름을 변수에 저장하고 싶습니다. 이름을 지정해 보겠습니다.
firstline
ㅏ두번째파일(Necha2_best_models.gff), 이 이름이 나타나는 모든 줄을 원합니다:
grep -ir "e_gw1.5.88.1" --include="Necha2*.gff" > Necha2_in_genome.list
그러나 명명된 변수의 경우:
grep -ir $firstline --include="Necha2*.gff" > Necha2_in_genome.list
이것은 "e_gw1.5.88.1"을 사용하여 작동합니다...
여기에서 생성한 파일은 내가 자르고 싶은 DNA 조각의 이름(sca_5_chr5_3_0)과 내가 원하는 부분(글리프 2390621에서 글리프 2392655까지)을 알려줍니다. 세 번째 파일에서 올바른 비트를 얻으려면 이 정보가 필요합니다.
a=(`awk -F '\t' '{print $4}' Necha2_in_genome.list`)
startDNA=${a[1]}
endDNA=${a[${#a[@]}-1]}
#add 1000 or other number, depending on the problems with the gene
correctedstartDNA=$(($startDNA-1000))
correctedendDNA=$(($endDNA-1000))
안에제삼이 예에서는 다음에서 키워드(sca_5_chr5_3_0) 뒤의 일부 부분을 잘라내고 싶습니다.
Kamaraj와 hschou 덕분에 이제 이에 대한 부분적인 해결책이 생겼습니다.
cat Necha2_scaffolds.fasta | sed -n -e '/sca_5_chr5_3_0/,$p' | grep -v '^>' | tr -d '\n'|awk -v start="${correctedstartDNA}" -v end="${correctedendDNA}" '{print substr($0,start,end)}' RS= Necha2_scaffolds.fasta
그러나 더 작은 숫자로 디버깅하는 경우:
cat Necha2_scaffolds.fasta | sed -n -e '/sca_5_chr5_3_0/,$p' | grep -v '^>' | tr -d '\n'|awk -v start=10 -v end=20 '{print substr($0,start,end)}' RS= Necha2_scaffolds.fasta
나는 다음과 같은 결과를 얻습니다.
r1_1_0
CCTTATCCTAGCG
nmapped
CTTATATATTAT
nmapped
TAAAAGGAGTTA
unmapped
TCTTATATAAA
unmapped
AATCTTAAGAA
RS 옵션이 무시되고 10~20의 특정 줄만 인쇄되는 것 같습니다. 왜 이 라인이 선택되었는지 모르겠습니다.
sca_5_chr5_3_0
파일에 한 번만 나타납니다.
다른 이름도 있어요
>sca_66_unmapped
>sca_67_unmapped
등
나는 178개의 게놈에서 이 정보를 얻어야 하는데, 그것들은 거대한 파일이고 수동 검색은 옵션이 아닙니다.
파일 모양은 다음과 같습니다.
N_haematococca_1.fasta(파일 1)는 일반 fasta 파일입니다.
>jgi|Necha2|40463|e_gw1.5.88.1
MTTRAPKDTYAQPDYEEAHLATFAAPKGYPIESMLPPSVKRETFEQALAEFTDAIGKDYVFIGDALSHYI
DPYDIYIDDSERRRMPSAAVCPSSLEEAKQALKVANKYGIPIWAFSRGKNLGYGGPSARVNGSVAFDLHR
Necha2_best_models.gff(파일 2)는 다음과 같습니다(더 길다):
sca_5_chr5_3_0 JGI exon 2390621 2390892 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2390621 2390892 . + 0 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 1
sca_5_chr5_3_0 JGI start_codon 2390621 2390623 . + 0 name "e_gw1.5.88.1"
sca_5_chr5_3_0 JGI exon 2390949 2391041 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2390949 2391041 . + 2 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 2
sca_5_chr5_3_0 JGI exon 2391104 2391311 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2391104 2391311 . + 2 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 3
sca_5_chr5_3_0 JGI exon 2391380 2392367 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2391380 2392367 . + 0 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 4
sca_5_chr5_3_0 JGI exon 2392421 2392485 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2392421 2392485 . + 1 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 5
sca_5_chr5_3_0 JGI exon 2392541 2392657 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2392541 2392657 . + 0 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 6
sca_5_chr5_3_0 JGI stop_codon 2392655 2392657 . + 0 name "e_gw1.5.88.1"
sca_5_chr5_3_0 JGI exon 2396205 2396730 . + . name "e_gw1.5.227.1"; transcriptId 41333
sca_5_chr5_3_0 JGI CDS 2396205 2396730 . + 0 name "e_gw1.5.227.1"; proteinId 41333; exonNumber 1
sca_5_chr5_3_0 JGI start_codon 2396205 2396207 . + 0 name "e_gw1.5.227.1"
Necha2_scaffolds.fasta(파일 3)는 다음과 같습니다(더 긴 GATC 비트...).
>sca_8_chr1_1_0
CCTTATCCTAGCGAGGATTAAGGTCCTCGAAAAGAAAGGCAAAGATATAAAGGTAATATAATAGAAGATT
AAGGTATTCTAAGTAAAGGTTATAAAGAAATAAAATAAGAAGAATATTTATAGGCTAAGAAAGACCCCCC
TAAAGGTTAAGGACTTAATATTAAGATTTAATATTCCCTAATTAATTAATATATTAATAAAAATAAAGAT
>sca_5_chr5_3_0
ATGACCACTATATCCATCGGCACAACGGCGTTAATCACATTTGGGTCTGCAATTTTCTGTTTTTGCGGTT
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAGACTCGCTCTGCTTT
>sca_67_unmapped
CTACAGGAAGCCCTAGAGGCCCTGCAAGACCTTCCAAAGCAGCACTCTTTGCTTCTTCTTAGAGGCAGTA
TCCAGCTTCTTCTAAGGCACCTCCAGAGGCAGCTAGACCCAACCGGGCTAAAAGACCTTTGGGAAGAGGC
TAATACCCTTATAAGAGAGGCTATTATAGCCCTAGTGGCTAGAAGCCCTAGTGAGAGCCCTAAAGAGCCT
AATTCAAGCCTTATAGCCCTTCCAGTCAGAGAGGGAGGCCTAGGAATACCCTTACACAAGGACCTAGCCC
예상 최종 출력: >sca_5_chr5_3_0 키워드 뒤의 1비트 텍스트입니다.
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAG