Linux 生信項目案例在線分析


本文轉自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


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM