편집 및 솔루션
내 원래 질문의 표현이 잘못되어 바퀴를 다시 만들려고 했기 때문에 지금은 내 질문에 답하고 있습니다(다른 사람에게 도움이 될 수도 있음).
gff2fasta는 제가 필요한 작업을 정확하게 수행하는 도구입니다. 즉, 전체 게놈(FULLGENOME.fasta라는 fasta 형식의 거대한 파일)에서 특정 DNA 조각을 추출하는 것입니다.
내가 원하는 조각이 어디에 있는지 안다면, gff2fasta에 대한 내 조각의 비계(여기: sca_5_chr5_3_0)와 시작(여기: 2390621)과 끝(여기: 2391041)을 지정하는 TEST.gff라는 파일을 만들 수 있습니다. 예를 들어:
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과 다른 모든 분들께 감사드립니다!
-------------
이것은 더 많은 정보와 세부 정보가 포함된 이 질문의 연속입니다.유닉스: 파일에서 문자 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 키워드 뒤의 단일 텍스트 비트입니다.
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAG