Linux基礎(chǔ)(二)

緊接linux基礎(chǔ)(一)的內(nèi)容

11、安裝 miniconda,建立虛擬環(huán)境,在虛擬環(huán)境中安裝生信必備軟件

# miniconda官網(wǎng)下載linux版安裝包:Miniconda2-4.7.10-Linux-x86_64.sh*
>./Miniconda2-4.7.10-Linux-x86_64.sh*     #進入安裝包所在目錄,運行安裝包
>conda config --add channels r
>conda config --add channels conda-forge
>conda config --add channels bioconda
>conda create -n rna python=2  #建立python2環(huán)境的名叫rna的虛擬conda環(huán)境
>souce activate rna    #激活進入rna環(huán)境
>conda install -y bwa gatk4 sra-tools fastqc trim-galore star hisat2 bowtie2 subread htseq multiqc samtools salmon     #安裝軟件
>source deactivate  #退出環(huán)境

12、產(chǎn)生bam文件的命令,示例:

>bowtie2-align-s --wrapper basic-0 -p 20 -x /yanke/yanke21/database/index/bowtie2/hg38 -S WT_rep2_Input.sam -U /tmp/41440.unp
>samtools view tmp.sorted.bam | head -n 2     #運行虛擬環(huán)境中的軟件要先進入conda環(huán)境
SRR1042600.42157053     0       chr1    629895  42      51M     *       0       0       ATAACCAATACTACCAATCANTACTCATCATTAATAATCATAATGGCTATA      CCCFFFFFHHHHHJJJJJJJ#4AGHJJIIJJIIIIIJJJJIJIIIIJJIJI     AS:i:-6 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:11C8A30    YT:Z:UU
SRR1042600.42212881     0       chr1    629895  42      51M     *       0       0       ATAACCAATACTACCAATCANTACTCATCATTAATAATCATAATGGCTATA      @@<FDFFBFDHHFJEIIGJI#3AFHGEHEIJIIGIIGGIJIIJIGIIGIIJ     AS:i:-6 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:11C8A30    YT:Z:UU
SRR1042600.12010763     16      chr1    629895  24      51M     *       0       0       ATAACCAATACTTCTAATCAAAACTCATCATTAATAATCATAATGGCTATA      ?4B?1*4DD?11*1*?+22+<3F:3@EC:CC4EA,DEDDDDD?D3B:==+;     AS:i:-10        XN:i:0  XM:i:4  XO:i:0  XG:i:0  NM:i:4  MD:Z:11C0A1C6T29YT:Z:UU
SRR1042600.29629551     16      chr1    629895  40      51M     *       0       0       ATAACCAATACTACCAATCACTACTCATCATTAATAATCATAATGGCTATA      HGF?JJHHFDHHGJJIHDFA+E?JIJJIIHGJJJJJJJHHHHHFFFFFCC@     AS:i:-8 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:11C8A30    YT:Z:UU
SRR1042600.41910745     0       chr1    629896  42      51M     *       0       0       TAACCAATACTACCAATCAANACTCATCATTAATAATCATAATGGCTATAG      CC@FFFFFHHHHGIIHIJJJ#3<CFHCGGIIIJJJJJJJJIGGFHIIJFII     AS:i:-6 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:10C9T30    YT:Z:UU
SRR1042600.14329856     16      chr1    629896  8       18M1I32M        *       0       0       AAACCAAATCCTCCAATCAAATCCTCATCATTAATAATCATAATGGCTATA      #############################@IHHGCE9GHFHHHDDDDD<@@     AS:i:-18        XN:i:0  XM:i:5  XO:i:1  XG:i:1  NM:i:6  MD:Z:0T6T0A2A9A28        YT:Z:UU

13、查看使用的參考基因組 gencode.v29.annotation.gtf.gz 具體有多少條染色體

>zcat gencode.v29.annotation.gtf.gz | grep -v "#" | cut -f 1| sort | uniq -c
250524 chr1
  99843 chr10
 167540 chr11
 159530 chr12
  42767 chr13
  96124 chr14
 100426 chr15
 130878 chr16
 171137 chr17
  48848 chr18
 168227 chr19
 204338 chr2
  61485 chr20
  29652 chr21
  59764 chr22
 172809 chr3
 108328 chr4
 120515 chr5
 124117 chr6
 130671 chr7
  98583 chr8
  97204 chr9
    143 chrM
  89500 chrX
   9064 chrY

14、上面的后綴為BAM 的文件的第二列,只有 0 和 16 兩個數(shù)字,用 cut/sort/uniq等命令統(tǒng)計它們的個數(shù)。

>samtools view tmp.sorted.bam | cut -f 2 | sort | uniq -c 
   29 0
   24 16

15、重新打開 rmDuplicate/samtools/paired 文件夾下面的后綴為BAM 的文件,再次查看第二列,并且統(tǒng)計

>cd ~/rmDuplicate/samtools/paired
>ll
total 52
drwxrwxr-x 2 yanke21 yanke21  4096 Nov 12  2016 ./
drwxrwxr-x 4 yanke21 yanke21  4096 Nov 12  2016 ../
-rw-rw-r-- 1 yanke21 yanke21   354 Nov 12  2016 readme.txt
-rw-rw-r-- 1 yanke21 yanke21  3111 Nov 12  2016 tmp.header
-rw-rw-r-- 1 yanke21 yanke21  4379 Nov 12  2016 tmp.rmdup.bam
-rw-rw-r-- 1 yanke21 yanke21  1886 Nov 12  2016 tmp.rmdup.vcf.gz
-rw-rw-r-- 1 yanke21 yanke21 11649 Nov 12  2016 tmp.sam
-rw-rw-r-- 1 yanke21 yanke21  4543 Nov 12  2016 tmp.sorted.bam
-rw-rw-r-- 1 yanke21 yanke21  1888 Nov 12  2016 tmp.sorted.vcf.gz
>samtools view tmp.sorted.bam | cut -f 2 | sort | uniq -c 
     8 147
     3 163
     1 323
     1 353
     1 371
     1 387
     1 433
     3 83
     2 97
     9 99

16、下載 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件,并且解壓,查看里面的文件夾結(jié)構(gòu), 這個文件有2.3M,注意留心下載時間及下載速度。

>wget -c http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 
>unzip sickle-results.zip
>tree sickle-results
sickle-results
├── command.txt
├── single_tmp_fastqc.html
├── single_tmp_fastqc.zip
├── test1_fastqc.html
├── test1_fastqc.zip
├── test2_fastqc.html
├── test2_fastqc.zip
├── trimmed_output_file1_fastqc.html
├── trimmed_output_file1_fastqc.zip
├── trimmed_output_file2_fastqc.html
└── trimmed_output_file2_fastqc.zip

17、解壓 sickle-results/single_tmp_fastqc.zip 文件,并且進入解壓后的文件夾,找到 fastqc_data.txt 文件,并且搜索該文本文件以 >>開頭的有多少行?

>unzip sickle-results/single_tmp_fastqc.zip
>cd single_tmp_fastqc/ && ll
total 296
drwxrwxr-x  4 yanke21 yanke21   4096 Oct  6  2016 ./
drwxr-xr-x 13 yanke21 yanke     4096 Aug 24 11:29 ../
-rw-rw-r--  1 yanke21 yanke21  12048 Oct  6  2016 fastqc_data.txt
-rw-rw-r--  1 yanke21 yanke21   5186 Oct  6  2016 fastqc.fo
-rw-rw-r--  1 yanke21 yanke21 259290 Oct  6  2016 fastqc_report.html
drwxrwxr-x  2 yanke21 yanke21   4096 Oct  6  2016 Icons/
drwxrwxr-x  2 yanke21 yanke21   4096 Oct  6  2016 Images/
-rw-rw-r--  1 yanke21 yanke21    542 Oct  6  2016 summary.txt
>grep -E '^>>' fastqc_data.txt | wc -l
24

18、下載 http://www.biotrainee.com/jmzeng/tmp/hg38.tss文件,去NCBI找到TP53/BRCA1等自己感興趣的基因?qū)?yīng)的 refseq數(shù)據(jù)庫ID,然后找到它們的hg38.tss文件的哪一行。

https://www.ncbi.nlm.nih.gov/gene/7157

>wget -c  http://www.biotrainee.com/jmzeng/tmp/hg38.tss
>less hg38.tss
NM_012279       chr5    177020680       177024680       0
NR_131773       chr5    177020680       177024680       0
NR_131774       chr5    177020680       177024680       0
NM_001308213    chr5    177020680       177024680       0
NM_001308219    chr5    177020680       177024680       0
NM_001318082    chr7    99502644        99506644        0
NM_001291980    chr5    177084872       177088872       0
NM_001283039    chr11   3795819 3799819 0
NM_001283040    chr11   3795819 3799819 0
NM_138964       chr2    68641588        68645588        0
NM_139057       chr15   100340005       100344005       1
>grep 'NM_001143990' hg38.tss
NM_001143990    chr17   7684071 7688071 0
>grep 'NM_001291588' hg38.tss
NM_001291588    chr17   43209827        43213827        0

19、解析hg38.tss 文件,統(tǒng)計每條染色體的基因個數(shù)。

cut -f 2 hg38.tss | sort | uniq -c

20、解析hg38.tss 文件,統(tǒng)計NM和NR開頭的熟練,了解NM和NR開頭的含義

grep  -E '^NM' hg38.tss | wc -l
51064
grep  -E '^NR' hg38.tss | wc -l
15954

NM開頭表示轉(zhuǎn)錄產(chǎn)物序列,成熟的mRNA序列
NR開頭表示非編碼的轉(zhuǎn)錄子序列,包括RNAs,假基因轉(zhuǎn)子等

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容