本文轉自https://freeaihub.com/article/bioinformatics-in-linux.html,該頁可在線進行操作。
Linux對生物信息的學習和實踐有強大的輔助作用,不管自己編寫shell命令腳本處理數據,還是使用現有豐富生物信息分析工具。Linux強大命令行功能,可以快速、批量、靈活的處理數據的提取、統計和整理等耗時耗力的重復性工作。日常生信分析中,多數整理工作都是用Linux命令的組合完成的,這相比於寫完整的Python或Perl程序更簡便快捷。
本節我們將通過一個案例來了解如何使用Linux來進行Linux生項項目的實踐分析。
bowtie2簡介與在線安裝
Bowtie2 是將測序reads與長參考序列比對工具 (適用於將長度大約為50到100或1000字符的reads與相對較長的基因組, 如哺乳動物,進行比對)。
通常是比較基因組學(包括識別變體(variation calling),ChIP-seq,RNA-seq,BS-seq)管道的第一步。
可以處理非常長的讀數(即10s或100s的千字節),但它針對近期測序儀產生的讀數長度和誤差模式進行了優化,如Illumina HiSeq 2000,Roche 454和Ion Torrent儀器。
Bowtie2使用FM索引(基於Burrows-Wheeler Transform 或 BWT)對基因組進行索引,以此來保持其占用較小內存。
#安裝bowtie2
cd /usr/local
cp /share/tar/bowtie2-2.3.4.1-linux-x86_64.zip .
unzip bowtie2-2.3.4.1-linux-x86_64.zip
mv bowtie2-2.3.4.1-linux-x86_64 bowtie2
echo 'export PATH=/usr/local/bowtie2:$PATH' >> ~/.bashrc
source ~/.bashrc #生效環境變量
#驗證安裝是否成功
bowtie2 --version
samtool簡介與在線安裝
samtools是一個用於操作sam和bam文件的工具合集。能夠實現二進制查看、格式轉換、排序及合並等功能,結合sam格式中的flag、tag等信息,還可以完成比對結果的統計匯總。同時利用linux中的grep、awk等操作命令,還可以大大擴展samtools的使用范圍與功能。包含有許多命令。
#安裝samtool
cd /usr/local
cp /share/tar/samtools-1.9.tar.bz2 .
#wget -c https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar xvf samtools-1.9.tar.bz2
/usr/local/samtools-1.9/configure --prefix=/usr/local/samtools-1.9
echo 'export PATH=/usr/local/samtools-1.9:$PATH' >> ~/.bashrc
source ~/.bashrc #生效環境變量
samtools --version
下載相應分析數據
下載http://www.biotrainee.com/jmzeng/igv/test.bed 文件,后在里面選擇含有 H3K4me3 的那一行是第幾行,該文件總共有幾行。
pwd
wget –c http://www.biotrainee.com/jmzeng/igv/test.bed
ls
#(-n 標記行數,-o 只顯示匹配上的,--color匹配文字出現顏色)
grep -n -o --color H3K4me3 test.bed
#(wc顯示文件的行數、單詞數、字節數)
cat test.bed |wc -l
下載 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件,並且解壓,查看里面的文件夾結構
wget -c http://www.biotrainee.com/jmzeng/rmDuplicate.zip
ls
unzip rmDuplicate.zip
ls
tree rmDuplicate
打開解壓的文件,進入 rmDuplicate/samtools/single 文件夾里面,查看后綴為 .sam 的文件,搞清楚生物信息學里面的SAM/BAM定義是什么【另外關於sam頭部注釋信息在tmp.header中】
SAM的全稱是sequence alignment/map format, 主要用於測序序列mapping到基因組上的結果表示
BAM是SAM的二進制壓縮文件,不能像sam可以用less查看,它需要用samtools view打開
SAM分數據通常有兩部分,頭部注釋信息(header section )和主體比對結果部分 (alignment section)
NM:i 經過編輯的序列長度(edit distance)
MD:Z 錯配位置/鹼基(mismatching positions/bases)
AS:i 匹配得分(Alignment score)
XS:i 第二好的匹配得分(suboptimal alignment score)
YS:i mate序列匹配的得分
BC: 條碼序列(barcode sequence)
rmDuplicate/samtools/single中的sorted.bam文件中,統計第二列各flag個數
cd ~/rmDuplicate/samtools/single
samtools view tmp.sorted.bam | cut -f2 | sort -n | uniq -c
下載 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件,並且解壓,進入解壓后的文件夾test1_fastqc,找到 fastqc_data.txt 文件,統計 以>>開頭的有多少行
cd ~
wget http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
unzip sickle-results.zip
cd sickle-results
unzip test1_fastqc.zip
cd test1_fastqc
cat fastqc_data.txt | grep '>>' | wc -l
下載 http://www.biotrainee.com/jmzeng/tmp/hg38.tss (存儲轉錄起始位點信息)文件,在NCBI中找到自己感興趣的基因對應的 refseq數據庫 ID(這里我以乳腺癌BRCA1為例),然后找到它在hg38.tss 文件的哪一行
cd ~
wget http://www.biotrainee.com/jmzeng/tmp/hg38.tss
cat hg38.tss | grep 'NM_007294'
關於RefSeq數據庫:
RefSeq數據庫中所有的數據是一個非冗余的、提供參考標准的數據,包括染色體、基因組(細胞器、病毒、質粒)、蛋白、RNA等。refseq序列是NCBI篩選過的非冗余數據庫,比GeneBank准確。
關於它的ID:NM開頭的表示標准序列,XM表示預測的蛋白編碼序列,NR表示非編碼蛋白的mRNA序列,AF開頭的表示克隆序列,BC開頭的表示模板序列
另外,你可能見過gi|4557284|ref|NM_000646.1|[4557284]這種格式
gi就是代表genebank identifier;ref就是對應的refseq中的ID啦
解析hg38.tss 文件,統計每條染色體的基因個數
cd ~
cut -f2 hg38.tss |cut -d'_' -f1 | sort |uniq -c | sort -rn
統計hg38.tss 文件中NM和NR開頭的序列,了解NM和NR開頭的含義
#統計NM開頭
cat hg38.tss | grep 'NM' | wc -l
#同理可以統計NR開頭
#NM開頭的表示標准序列,可以轉錄成蛋白質的基因
#NR非編碼蛋白的mRNA序列
本節總結了兩款生信專用軟件bowtie2和samtool的安裝,並使用了Linux中的命令,如wc
,cut
等常見命令完成生物信息項目的簡單分析。
參考:1.劉小澤 《Linux學習》 簡書 https://www.jianshu.com/p/0b0e7c0f31db