【생물정보학】 리눅스 프로그래밍(bash programming)

리눅스 프로그래밍(bash programming)

 

추천글 : 【운영체제】 리눅스(Linux) 

---

1. 문법 [본문]

2. 생물정보학 예제 [본문]

---

 

1. 문법 [목차]

⑴ 변수(variable)

① 개요 : 긴 string을 나타내거나 loop와 같은 bash programming에 유용함

② 변수 지정 : $ 기호를 사용 

 

i="Hello World!"
echo $i

 

③ 모호함을 줄이기 위해 중괄호를 사용할 수 있음

 

i="a"
i¡="b"
echo $ii (will give b)
echo ${i}i (will give ai)

 

⑵ for 루프

 

for varname in list:
do 
    command1
    command2
    ..
done

 

 

2. 생물정보학 예시 [목차]

⑴ wc -l filename 

① 해석 : 라인 수 세기

② wildcard(wc)를 rm, cp, mv와 함께 쓸 때는 유의해야 함

⑵ grep -v ^# gencode.v19.long_noncoding_RNAs.gtf | wc -l 

① ^ (caret)의 역할 : ^는 정규표현식에서 줄의 시작을 나타냄. ^#는 "줄의 시작이 #로 시작하는 줄"을 의미

② -v의 역할 :-v는 반대를 뜻함. -v를 쓰지 않으면, 반대로 #로 시작하는 줄만 출력됨

③ | (vertical bar)의 역할 : 앞의 명령어의 출력 결과를 다음 명령어인 wc - l의 입력으로 전달

④ grep -v ^# gencode.v19.long_noncoding_RNAs.gtf : 파일에서 #로 시작하는 줄을 제외한 나머지 줄을 모두 검색

⑤ wc -l  : 이전 커멘드에서 넘겨준 줄의 개수를 출력

⑥ 해석 : 파일에서 #로 시작하는 줄을 제외한 나머지 줄의 개수를 출력

⑶ grep -w gene gencode.v19.long_noncoding_RNAs.gtf | wc -l 

① 해석 : 파일에서 독립적인 단어 gene이 포함된 줄의 개수를 출력

② -w 없이 grep gene을 사용하면, gene이 부분적으로 포함된 모든 줄이 검색 : gene, genetics, progene 등

⑶ for i in `seq 1 3`; do mkdir sample${i}; done

① backtick(`)의 역할 : ` ` 안에 있는 커멘드를 실행하고 얻은 아웃풋을 나타냄

⑷ grep "^chr6" gencode.v19.long_noncoding_RNAs.gtf | perl -lane 'if ($F[3]>=28000000 & ($F[4]<=34000000)){print $_;}' | grep -w "gene" 

① -l : 각 라인 끝에 있는 newline 문자를 제거함 

② -a : autosplit 모드. 즉, 띄어쓰기로 구분된 각 입력 라인을 @F array에 넣음

○ 첫 번째, 두 번째, ··· 원소를 각각 $F[0], $F[1], ···로 접근할 수 있음

○ tab으로 구분된 파일이라면, perl -F'\t' -lane 'print $F[1]' input.csv 와 같이 커멘드를 구성

③ -n : 각 라인이 $_으로 나타내어짐 

④ -e : perl 프로그램이 커멘드 라인에 들어갈 수 있게 함

⑤ grep "^chr6" gencode.v19.long_noncoding_RNAs.gtf : chr6로 시작하는 줄을 모두 검색

⑥ perl -lane 'if ($F[3]>=28000000 & ($F[4]<=34000000)){print $_;}' : perl 프로그래밍을 커멘드 라인으로 구현하고, 3번째 값 (e.g., start)이 28M 이상이고 4번째 값 (e.g., end)이 34M 이하인 라인을 모두 검색

⑦ grep -w "gene" : 이전 커멘드에서 넘겨 준 라인에서 "gene"에 해당하는 부분을 검색  

⑧ 해석 : chr6 28-34M에 있는 유전자 탐색

⑸ grep -w "gene" genecode.v19.long_noncoding_RNAs.gtf | perl -lane 'm/gene_name \"(.*?)\"/; print $1' | less -S 

① m/gene_name \"(.*)\"/ : 패턴 매칭

② m/gene_name \"(.*)\"/ : 따옴표 처리

③ m/gene_name \"(.*)\"/ : capture group으로 $1로 호출할 수 있음

○ . : 줄바꿈 문자를 제외한 아무 문자와 매칭됨. a, 1, # 등 어떤 문자와도 매칭됨

○ * : 바로 앞의 문자나 패턴이 0번 이상 반복될 수 있음을 나타냄

○ a* : 빈 문자열(""), a, aa, aaa 등과 매칭됨 

○ .* : 아무 문자열과 매칭됨 (0번 이상 반복)

○ "(.*?)" : non-greedy match. 즉, 큰따옴표 안에 있는 텍스트를 최소한의 문자만 포함하도록 매칭. ?이 없으면 라인 끝까지 검색

⑹ j='echo $i | sed 's/.bam//g"

① sed : '  '를 통해 정규표현식을 나타냄

② s/.bam// : 치환(substitution). s/패턴/변경할_문자열/변경될_문자열 

③ s/.bam// : 변경될 문자열이 빈 문자열이므로 .bam이 삭제됨 

④ g : 한 줄에서 .bam이라는 문자열이 여러 번 나타날 경우, 모든 .bam을 제거. 만약 g 플래그가 없으면, 가장 첫 번째 .bam만 제거

⑺ samtools view -h $i | perl -lane 'if (($F[0]=~/^@/ | ($F[6] eq"=")) {print $_;}' | samtools view -S -b - > $(j}.new.bam

① samtools 

○ -h : SAM 헤더를 포함하도록 지정 

○ -S : 입력이 SAM 포맷임을 의미함. samtools 최신 버전에서는 이를 요하지 않음

○ -b : 출력을 BAM 포맷으로 지정 

② $F[0]=~ /^@/ : $F[0]이 @로 시작하는 경우(헤더 라인)

③ $F[6]이 "="인 경우(자기 자신에게 정렬된 라인)

④ 해석 : paired mapping에서 다른 염색체에 맵핑된 리드를 제거

 

입력: 2024.01.16 00:28