原本應(yīng)該是學(xué)習(xí)的時(shí)候記錄的，但是各種原因耽誤了，只能想起來了就記一點(diǎn)

比對

可以使?用bowtie2軟件?自帶的測試數(shù)據(jù)?生成sam/bam?文件，代碼如下:

mkdir -p ~/biosoft
cd ~/biosoft
# 針對不不同系統(tǒng)下載不不同的版本
# 優(yōu)先下載最新版
wget https://sourceforge.net/projects/bowtie- bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip
unzip bowtie2-2.3.4.3-linux-x86_64.zip
cd ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads
../../bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq > tmp.sam # samtools view -bS tmp.sam >tmp.bam

也可以用conda 直接安裝bowtie2和samtools軟件。

(base) ChengLiangPing@VM-0-8-ubuntu:~$ ls
11-17  Miniconda2-latest-Linux-x86_64.sh  biosoft  miniconda2  readme.txt
(base) ChengLiangPing@VM-0-8-ubuntu:~$ cd ~/biosoft/
(base) ChengLiangPing@VM-0-8-ubuntu:~/biosoft$ ls
Miniconda2-latest-Linux-x86_64.sh
(base) ChengLiangPing@VM-0-8-ubuntu:~/biosoft$ conda install -y bowtie2 samtoolsCollecting package metadata (current_repodata.json): done
Solving environment: done

找變異位點(diǎn)

mkdir -p ~/biosoft
cd ~/biosoft
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-
latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.bashrc
conda config --add channels
https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels
https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels
https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
conda  create -y -n test
conda activate test
conda install -y samtools bcftools
wget https://sourceforge.net/projects/bowtie-
bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip
unzip bowtie2-2.3.4.3-linux-x86_64.zip
cd ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads
../../bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq | samtools
sort -@ 5 -o tmp.bam -
bcftools mpileup -f ../reference/lambda_virus.fa tmp.bam |bcftools call -vm >
tmp.vcf

然后把這些?文件，全部載?入IGV 進(jìn)?行行可視化。

作業(yè):完成?一個(gè)SNP-calling模擬項(xiàng)?

首先下載X,Y染?色體的fasta序列列，在UCSC上?面下載即可。
然后把X染?色體構(gòu)建bwa的索引 接著模擬?一個(gè)Y染?色體的測序數(shù)據(jù)，模擬的程序很簡單,模擬Y染?色體的測序?片段(PE100， insert400)
然后把模擬測序數(shù)據(jù)?比對到X染?色體的參考，統(tǒng)計(jì)?一下比對結(jié)果。 最后對比對成功的bam?件進(jìn)行找變異位點(diǎn)。
代碼如下:

## 源代碼?方式安裝 bwa-0.7.15
## conda安裝samtools
cd tmp/chrX_Y/hg19/
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrX.fa.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrY.fa.gz
gunzip chrX.fa.gz
gunzip chrY.fa.gz
 
~/biosoft/bwa/bwa-0.7.15/bwa index chrX.fa
wget https://raw.githubusercontent.com/jmzeng1314/my-perl/master/2.chrX- chrY/simulate.pl
perl simulate.pl chrY.fa ## 這個(gè)perl腳本并不不需要大家理解。
~/biosoft/bwa/bwa-0.7.15/bwa mem -t 5 -M chrX.fa read*.fa >read.sam
samtools view -bS read.sam >read.bam
samtools flagstat read.bam
samtools sort -@ 5 -o read.sorted.bam read.bam
samtools view -h -F4 -q 5 read.sorted.bam |samtools view -bS|samtools rmdup - read.filter.rmdup.bam
samtools index read.filter.rmdup.bam
# 后面bam-->vcf代碼?己嘗試。
# samtools mpileup -ugf ~/tmp/chrX_Y/hg19/chrX.fa read.filter.rmdup.bam |bcftools call -vmO z -o read.bcftools.vcf.gz
## 把fa/bam/vcf 載?入到 IGV 進(jìn)?行行可視化，截圖其中?一個(gè)變異位點(diǎn)
## 參考 http://www.biotrainee.com/thread-696-1-1.html

一個(gè)RNA項(xiàng)?分析流程

通常，RNA項(xiàng)?指的是轉(zhuǎn)錄組測序，其研究對象為特定細(xì)胞群體在某一功能狀態(tài)下所能轉(zhuǎn)錄出來的所有 RNA 的總和，包括 mRNA 和非編碼 RNA 。通過轉(zhuǎn)錄組測序，能夠全?面獲得物種特定組織或器官的轉(zhuǎn)錄本信息，從?進(jìn)行轉(zhuǎn)錄本結(jié)構(gòu)研究、變異研究、基因表達(dá)?平研究以及全新轉(zhuǎn)錄本發(fā)現(xiàn)等研究。

其中，基因表達(dá)?平的探究是轉(zhuǎn)錄組領(lǐng)域最熱?的方向，利用轉(zhuǎn)錄組數(shù)據(jù)來識別轉(zhuǎn)錄本和表達(dá)定量，是轉(zhuǎn)錄組數(shù)據(jù)的核?作?。由于這個(gè)作用，他可以不不依賴其他組學(xué)信息，單獨(dú)成為一個(gè)產(chǎn)品項(xiàng)?：RNA-seq 測序。所以很多時(shí)候轉(zhuǎn)錄組測序會與RNA-seq混為一談。
現(xiàn)在RNA-seq數(shù)據(jù)使?廣泛，但是沒有一套流程可以解決所有的問題。

RNA-seq分析中的重要的步驟包括:實(shí)驗(yàn)設(shè)計(jì)，質(zhì)控，read比對，表達(dá)定量，可視化，差異表達(dá)，識別可變剪切，功能注釋，融合基因檢測，eQTL定位等。

來自于 R處理 mRNA-seq數(shù)據(jù)教程:
http://biocluster.ucr.edu/~rkaundal/workshops/R_mar2016/RNAseq.html

關(guān)于數(shù)據(jù)格式的信息，可以參考生信技能樹的數(shù)據(jù)格式專題，也可以參考UCSC的介紹：http://genome.ucsc.edu/FAQ/FAQformat.html

建議：請放棄GTF文件格式，改用GenePred格式處理轉(zhuǎn)錄組注釋！ http://www.biotrainee.com/thread-928-1-1.html

理解bw,bigwig文件是簡化的bam文件，http://www.biotrainee.com/thread-797-1-1.html 參考這個(gè)【直播】我的基因組 35：bam格式轉(zhuǎn)化為bw格式看測序深度分布。

bedtools用法大全，它可以替代比較差的生信工程師：
http://www.biotrainee.com/thread-153-1-1.html

• shell中的擴(kuò)展(Expansions) https://opengers.github.io/linux/linux-shell-brace-parameter-command-pathname-expansion/
• bash腳本的參數(shù)擴(kuò)展 (parameter expansion) ：https://www.ibm.com/developerworks/cn/linux/l-bash-parameters.html
• shell通配符（wildcard）： https://cloud.tencent.com/developer/article/1114732
• type命令：https://man.linuxde.net/type
• 字符串操作：https://my.oschina.net/aiguozhe/blog/41557

shell擴(kuò)展有以下幾種，并按以下順序處理，當(dāng)然如果沒找到匹配的擴(kuò)展格式，那就不處理

brace expansion 花括號({})擴(kuò)展
tilde expansion ~字符擴(kuò)展
parameter and variable expansion 參數(shù)和變量擴(kuò)展
arithmetic expansion 算術(shù)擴(kuò)展
command substitution 命令替換
process substitution 過程替換
word splitting
Filename Expansion 通配符擴(kuò)展

解讀長指令：

less -S example.gtf  | awk -v FS="\t" '$3 == "gene"{print $9}' |head -n 5
less -S example.gtf  | awk -v FS="\t" '$3 == "gene"{split($9, x,";");print x[1]}' |head -n 5
less -S example.gtf  | awk -v FS="\t" '$3 == "gene"{split($9, x,";");split(x[1],y," ");print y[2]}' |head -n 5
less -S example.gtf  | awk -v FS="\t" '$3 == "gene"{split($9, x,";");split(x[1],y," ");gsub("\"","",name);print name}' |head -n 5
less -S example.gtf  | awk -v FS="\t" '$3 == "gene"{split($9, x,";");split(x[1],y," ");name=y[2];gsub("\"","",name);print name}' |head -n 5

使用hisat2來?對RNA-seq數(shù)據(jù)

hisat2和bowtie2都是約翰霍普金斯醫(yī)學(xué)院的科研?員開發(fā)的，所以用法類似。
?如公眾號推的: 玩轉(zhuǎn)RNA-seq數(shù)據(jù)也可以不不需要linux ?https://mp.weixin.qq.com/s/HzLbrBuEwR0cKwPZj-1Fkw 就介紹過一個(gè)數(shù)據(jù)集:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE109480

下載配套參考基因組的gtf?文件

?或者小鼠，推薦使用gencode數(shù)據(jù)庫 https://www.gencodegenes.org/ ，其它常見物種，推薦使用 ENSEMBL :https://asia.ensembl.org/info/data/ftp/index.html

使用featurecounts定量

根據(jù)正確的gtf?件對我們前面步驟的比對得到的bam文件進(jìn)行計(jì)數(shù)!
常?的基于?對的基因定量軟件:Htseq-count，bedtools mutilcov，featureCount，因?yàn)閭€(gè)人習(xí)慣， 講師這?著重講解 featureCount的用法。
featureCount是subread套件的?個(gè)模塊，最?的優(yōu)點(diǎn)就是速度?常快，使用全部overlap的reads計(jì) 數(shù)，靈活考慮多比對的reads的計(jì)數(shù)。 subread 軟件官?網(wǎng):http://subread.sourceforge.net/ http://bioi nf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf

featurecounts軟件的主要參數(shù):

-a 輸?入GTF/GFF基因組注釋?文件
-p 這個(gè)參數(shù)是針對paired-end數(shù)據(jù)
-F 指定-a注釋?文件的格式，默認(rèn)是GTF
-g 從注釋?文件中提取Meta-features信息?用于read count，默認(rèn)是gene_id -t 跟-g?一樣的意思， 其是默認(rèn)將exon作為?一個(gè)feature
-o 輸出?文件
-T 多線程數(shù)

featurecounts軟件運(yùn)行通常設(shè)置: -T 5 -p -t exon -g gene_id -a hg19.gtf -o out.txt input.bam

GTF的產(chǎn)生和流行有其歷史的原因。但是從技術(shù)角度來講，這個(gè)文件格式是個(gè)非常差勁的格式。

GTF格式非常冗余。以人類轉(zhuǎn)錄組為例，Gencode V22的GTF文件為1.2G，壓縮之后只有40M。大家知道壓縮軟件的壓縮比是和軟件的冗余程度。很少有文件能夠壓縮到1/30的大小?？梢奊TF格式多么冗余。GTF格式（及其早期版本GFF等）有很好的替代格式。從信息量上來講：GTF 等價(jià)于 GenePred （Bed12) + Gene_Anno_table。GenePred是Jimmy Kent創(chuàng)建UCSC genome browser的時(shí)候建立的文件格式。UCSC的文件格式定義是非常smart的，包括之后我可能會講到的2bit，bigwig格式。

GTF vs GenePred：
從文件大小上來看，壓縮前：GTF（1.2G) >> Genepred(23M) + Gene_Anno_table (2.8M)。壓縮后：GTF(40M) >> GenePred(7.8M) +Gene_Anno_table (580K)
從可讀性來講，GTF是以gene interval 為單位（行），每行可以是gene，transcript，exon，intron，utr等各種信息，看起來什么都在里面，很全面，其實(shí)可讀性非常差，而且容易產(chǎn)生各種錯(cuò)誤。GenePred格式是以transcript為單位，每行就是一個(gè)transcript，簡潔直觀。
從程序處理的角度來講：以GTF文件作為輸入的程序，如果換成以GenePred格式為輸入，編程的難度會降低一個(gè)數(shù)量級，運(yùn)算時(shí)間會快很多，代碼的可讀性強(qiáng)很多。

GTF 轉(zhuǎn)換成GenePred：

gtfToGenePred -genePredExt -ignoreGroupsWithoutExons -geneNameAsName2 test.gtf test.gpd

Gene_Anno_table: 其實(shí)就是把GTF的所有transcript行的第9列轉(zhuǎn)換變成一個(gè)表格。

應(yīng)用實(shí)例：
生信編程直播第一題：人類基因組的外顯子區(qū)域到底有多長？
1. 取出所有g(shù)ene的exon。
2. 對exon進(jìn)行排序。
3. 對有overlap的exon進(jìn)行merge。
4. 計(jì)算merge后的exon長度。
代碼如下：

def cmpBed(x,y):
    return cmp(x[0],y[0]) or cmp(x[1],y[1])
def isOverlap(bed1,bed2):
    return not (bed1[2]!=bed2[2] or bed1[0] > bed2[1] or bed1[1]<bed2[0])
def mergeBed(beds):
    tbed = beds[0]
    for bed in beds[1:]:
        if isOverlap(bed,tbed):
            tbed = (min(bed[0],tbed[0]),max(bed[1],tbed[1]),tbed[2])
        else:
            yield tbed
            tbed = bed
    yield tbed
 
exons = []
for line in open("test.genepred"):
    gene = line.rstrip().split('\t')
    exons.extend([(int(s),int(e),gene[1]) for s,e in zip(gene[8].rstrip(',').split(','),gene[9].rstrip(',').split(','))])
exons = sorted(exons,cmp=cmpBed)
print sum([bed[1]-bed[0] for bed in mergeBed(exons)])

生信編程直播第三題：hg38每條染色體基因,轉(zhuǎn)錄本的分布
1. 讀取genepred格式文件為DataFrame。
2. 按照chrom進(jìn)行g(shù)roup，然后count，最后barplot。
3. 按照gene symbol去重復(fù)，然后按照chrom進(jìn)行g(shù)roup，然后count，最后barplot。

import pandas
 
df = pandas.read_table("/scratch/bcb/ywang52/TData/genomes/hg38/gencode.v22.annotation_GeneSymbol.gpd.gz",header=None)
# transcript count
df.ix[:,:1].groupby(1).count().plot.bar(legend=False)
ylabel('Transcript count')
title('Gencode V22')
 
# gene count
sdf = df.ix[:,[1,11]].drop_duplicates().groupby(1).count().plot.bar(legend=False)
ylabel('Gene count')
title('Gencode V22')

生信編程直播第四題：多個(gè)同樣的行列式文件合并起來

雖然這個(gè)題目與genepred無關(guān)，但是還是做了吧。
1. 把每個(gè)文件讀成dataframe，用第一列做行名, 文件名做列名。
2. 按照行名merge dataframe

import os
import pandas
dfs = []
for f in os.listdir("."):
    if f.endswith(".readcount"):
        df = pandas.read_table(f,index_col=0,header=None)
        df.columns = [os.path.splitext(f)[0]] # file name as column name
        dfs.append(df)
data = pandas.concat(dfs,axis=1) # merge dataframe by row names

生信編程直播第五題：根據(jù)GTF畫基因的多個(gè)轉(zhuǎn)錄本結(jié)構(gòu)
以ANXA1基因?yàn)槔?br> 1. 按行讀取genepred文件，第3，4列為轉(zhuǎn)錄本的區(qū)間，第4，5列為ORF的區(qū)間，第9和10列為exon起始和終止位置。
2. 先畫Intron，即基因全長：第3，4列。
3. 再畫ORF 和UTR。把第9，10列分割成exon 的starts，stops。遍歷每個(gè)exon：
如果exon和ORF沒有overlap：畫成UTR。
反之： 先畫成UTR，再把overlap的部分畫成ORF。
NOTE: 這里有個(gè)trick，后畫的部分會覆蓋先畫的部分。

import matplotlib.pyplot as pltcolors = ['red','blue','purple']  # color loops
gids = []
for cnt,line in enumerate(open("ANXA1.genepred")):
    gene = line.rstrip().split("\t")
    gids.append(gene[0])
    plt.plot([int(gene[3]),int(gene[4])],[cnt,cnt],linewidth=1,color='black') # draw intron
    txstart, txstop = int(gene[5]),int(gene[6])
    for start,stop in zip([int(s) for s in gene[8].rstrip(',').split(',')],[int(s) for s in gene[9].rstrip(',').split(',')]):
        if start > txstop or stop < txstart: # UTR: no overlap with ORF
            plt.plot([start,stop],[cnt,cnt],linewidth=3,color='grey')
        else: # draw ORF
            plt.plot([start,stop],[cnt,cnt],linewidth=3,color='grey')
            plt.plot([max(txstart,start),min(txstop,stop)],[cnt,cnt],linewidth=5,color=colors[cnt%len(colors)])                
plt.yticks(range(cnt+1),gids)
plt.ylim(-0.5,cnt+0.5)
plt.tight_layout()
plt.savefig("test.pdf")

以GenePred格式為輸入的程序，所有代碼都不到20行，而且程序清晰簡潔，可讀性強(qiáng)（當(dāng)然Python語言的支持也是功不可沒），對我等新手非常友好 。
附：
UCSC 格式匯總: https://genome.ucsc.edu/FAQ/FAQformat
gtfToGenePred 二進(jìn)制文件： http://hgdownload.cse.ucsc.edu/a ... 86_64/gtfToGenePred
以上內(nèi)容整理自生信技能樹。